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 double LinearAlgebra::calcKruskalWallis(vector<spearmanRank>& values, double& pValue){
 810         try {
 811         double H;
 812         set<string> treatments;
 813
 814         //rank values
 815         sort(values.begin(), values.end(), compareSpearman);
 816         vector<spearmanRank*> ties;
 817         int rankTotal = 0;
 818         vector<int> TIES;
 819         for (int j = 0; j < values.size(); j++) {
 820             treatments.insert(values[j].name);
 821             rankTotal += (j+1);
 822             ties.push_back(&(values[j]));
 823
 824             if (j != values.size()-1) { // you are not the last so you can look ahead
 825                 if (values[j].score != values[j+1].score) { // you are done with ties, rank them and continue
 826                     if (ties.size() > 1) { TIES.push_back(ties.size()); }
 827                     for (int k = 0; k < ties.size(); k++) {
 828                         double thisrank = rankTotal / (double) ties.size();
 829                         (*ties[k]).score = thisrank;
 830                     }
 831                     ties.clear();
 832                     rankTotal = 0;
 833                 }
 834             }else { // you are the last one
 835                 if (ties.size() > 1) { TIES.push_back(ties.size()); }
 836                 for (int k = 0; k < ties.size(); k++) {
 837                     double thisrank = rankTotal / (double) ties.size();
 838                     (*ties[k]).score = thisrank;
 839                 }
 840             }
 841         }
 842
 843
 844         // H = 12/(N*(N+1)) * (sum Ti^2/n) - 3(N+1)
 845         map<string, double> sums;
 846         map<string, double> counts;
 847         for (set<string>::iterator it = treatments.begin(); it != treatments.end(); it++) { sums[*it] = 0.0; counts[*it] = 0; }
 848
 849         for (int j = 0; j < values.size(); j++) {
 850             sums[values[j].name] += values[j].score;
 851             counts[values[j].name]+= 1.0;
 852         }
 853
 854         double middleTerm = 0.0;
 855         for (set<string>::iterator it = treatments.begin(); it != treatments.end(); it++) {
 856             middleTerm += ((sums[*it]*sums[*it])/counts[*it]);
 857         }
 858
 859         double firstTerm = 12 / (double) (values.size()*(values.size()+1));
 860         double lastTerm = 3 * (values.size()+1);
 861
 862         H = firstTerm * middleTerm - lastTerm;
 863
 864         //adjust for ties
 865         if (TIES.size() != 0) {
 866             double sum = 0.0;
 867             for (int j = 0; j < TIES.size(); j++) { sum += ((TIES[j]*TIES[j]*TIES[j])-TIES[j]); }
 868             double result = 1.0 - (sum / (double) ((values.size()*values.size()*values.size())-values.size()));
 869             H /= result;
 870         }
 871
 872         //Numerical Recipes pg221
 873         pValue = 1.0 - (gammp(((treatments.size()-1)/(double)2.0), H/2.0));
 874
 875         return H;
 876     }
 877         catch(exception& e) {
 878                 m->errorOut(e, "LinearAlgebra", "calcKruskalWallis");
 879                 exit(1);
 880         }
 881 }
 882 /*********************************************************************************************************************************/
 883 //assumes both matrices are square and the same size
 884 double LinearAlgebra::calcKendall(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
 885         try {
 886                 double r;
 887
 888                 //format data
 889                 vector<spearmanRank> scores;
 890                 for (int i = 0; i < euclidDists.size(); i++) {
 891                         for (int j = 0; j < i; j++) {
 892                                 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
 893                                 scores.push_back(member);
 894                         }
 895                 }
 896
 897                 //sort scores
 898                 sort(scores.begin(), scores.end(), compareSpearman);
 899
 900                 //find ranks of xi
 901                 map<string, float> rankEuclid;
 902                 vector<spearmanRank> ties;
 903                 int rankTotal = 0;
 904                 for (int j = 0; j < scores.size(); j++) {
 905                         rankTotal += (j+1);
 906                         ties.push_back(scores[j]);
 907
 908                         if (j != (scores.size()-1)) { // you are not the last so you can look ahead
 909                                 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
 910
 911                                         for (int k = 0; k < ties.size(); k++) {
 912                                                 float thisrank = rankTotal / (float) ties.size();
 913                                                 rankEuclid[ties[k].name] = thisrank;
 914                                         }
 915                                         ties.clear();
 916                                         rankTotal = 0;
 917                                 }
 918                         }else { // you are the last one
 919
 920                                 for (int k = 0; k < ties.size(); k++) {
 921                                         float thisrank = rankTotal / (float) ties.size();
 922                                         rankEuclid[ties[k].name] = thisrank;
 923                                 }
 924                         }
 925                 }
 926
 927                 vector<spearmanRank> scoresUser;
 928                 for (int i = 0; i < userDists.size(); i++) {
 929                         for (int j = 0; j < i; j++) {
 930                                 spearmanRank member(toString(scoresUser.size()), userDists[i][j]);
 931                                 scoresUser.push_back(member);
 932                         }
 933                 }
 934
 935                 //sort scores
 936                 sort(scoresUser.begin(), scoresUser.end(), compareSpearman);
 937
 938                 //find ranks of yi
 939                 map<string, float> rankUser;
 940                 ties.clear();
 941                 rankTotal = 0;
 942                 for (int j = 0; j < scoresUser.size(); j++) {
 943                         rankTotal += (j+1);
 944                         ties.push_back(scoresUser[j]);
 945
 946                         if (j != (scoresUser.size()-1)) { // you are not the last so you can look ahead
 947                                 if (scoresUser[j].score != scoresUser[j+1].score) { // you are done with ties, rank them and continue
 948
 949                                         for (int k = 0; k < ties.size(); k++) {
 950                                                 float thisrank = rankTotal / (float) ties.size();
 951                                                 rankUser[ties[k].name] = thisrank;
 952                                         }
 953                                         ties.clear();
 954                                         rankTotal = 0;
 955                                 }
 956                         }else { // you are the last one
 957
 958                                 for (int k = 0; k < ties.size(); k++) {
 959                                         float thisrank = rankTotal / (float) ties.size();
 960                                         rankUser[ties[k].name] = thisrank;
 961                                 }
 962                         }
 963                 }
 964
 965                 int numCoor = 0;
 966                 int numDisCoor = 0;
 967
 968                 //order user ranks
 969                 vector<spearmanRank> user;
 970                 for (int l = 0; l < scores.size(); l++) {
 971                         spearmanRank member(scores[l].name, rankUser[scores[l].name]);
 972                         user.push_back(member);
 973                 }
 974
 975                 int count = 0;
 976                 for (int l = 0; l < scores.size(); l++) {
 977
 978                         int numWithHigherRank = 0;
 979                         int numWithLowerRank = 0;
 980                         float thisrank = user[l].score;
 981
 982                         for (int u = l+1; u < scores.size(); u++) {
 983                                 if (user[u].score > thisrank) { numWithHigherRank++; }
 984                                 else if (user[u].score < thisrank) { numWithLowerRank++; }
 985                                 count++;
 986                         }
 987
 988                         numCoor += numWithHigherRank;
 989                         numDisCoor += numWithLowerRank;
 990                 }
 991
 992                 r = (numCoor - numDisCoor) / (float) count;
 993
 994                 //divide by zero error
 995                 if (isnan(r) || isinf(r)) { r = 0.0; }
 996
 997                 return r;
 998
 999         }
1000         catch(exception& e) {
1001                 m->errorOut(e, "LinearAlgebra", "calcKendall");
1002                 exit(1);
1003         }
1004 }
1005 /*********************************************************************************************************************************/
1006 double LinearAlgebra::calcKendall(vector<double>& x, vector<double>& y, double& sig){
1007         try {
1008                 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1009
1010                 //format data
1011                 vector<spearmanRank> xscores;
1012                 for (int i = 0; i < x.size(); i++) {
1013                         spearmanRank member(toString(i), x[i]);
1014                         xscores.push_back(member);
1015                 }
1016
1017                 //sort xscores
1018                 sort(xscores.begin(), xscores.end(), compareSpearman);
1019
1020                 //convert scores to ranks of x
1021                 vector<spearmanRank*> ties;
1022                 int rankTotal = 0;
1023                 for (int j = 0; j < xscores.size(); j++) {
1024                         rankTotal += (j+1);
1025                         ties.push_back(&(xscores[j]));
1026
1027                         if (j != xscores.size()-1) { // you are not the last so you can look ahead
1028                                 if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
1029                                         for (int k = 0; k < ties.size(); k++) {
1030                                                 float thisrank = rankTotal / (float) ties.size();
1031                                                 (*ties[k]).score = thisrank;
1032                                         }
1033                                         ties.clear();
1034                                         rankTotal = 0;
1035                                 }
1036                         }else { // you are the last one
1037                                 for (int k = 0; k < ties.size(); k++) {
1038                                         float thisrank = rankTotal / (float) ties.size();
1039                                         (*ties[k]).score = thisrank;
1040                                 }
1041                         }
1042                 }
1043
1044
1045                 //format data
1046                 vector<spearmanRank> yscores;
1047                 for (int j = 0; j < y.size(); j++) {
1048                         spearmanRank member(toString(j), y[j]);
1049                         yscores.push_back(member);
1050                 }
1051
1052                 //sort yscores
1053                 sort(yscores.begin(), yscores.end(), compareSpearman);
1054
1055                 //convert to ranks
1056                 map<string, float> rank;
1057                 vector<spearmanRank> yties;
1058                 rankTotal = 0;
1059                 for (int j = 0; j < yscores.size(); j++) {
1060                         rankTotal += (j+1);
1061                         yties.push_back(yscores[j]);
1062
1063                         if (j != yscores.size()-1) { // you are not the last so you can look ahead
1064                                 if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
1065                                         for (int k = 0; k < yties.size(); k++) {
1066                                                 float thisrank = rankTotal / (float) yties.size();
1067                                                 rank[yties[k].name] = thisrank;
1068                                         }
1069                                         yties.clear();
1070                                         rankTotal = 0;
1071                                 }
1072                         }else { // you are the last one
1073                                 for (int k = 0; k < yties.size(); k++) {
1074                                         float thisrank = rankTotal / (float) yties.size();
1075                                         rank[yties[k].name] = thisrank;
1076                                 }
1077                         }
1078                 }
1079
1080
1081                 int numCoor = 0;
1082                 int numDisCoor = 0;
1083
1084                 //associate x and y
1085                 vector<spearmanRank> otus;
1086                 for (int l = 0; l < xscores.size(); l++) {
1087                         spearmanRank member(xscores[l].name, rank[xscores[l].name]);
1088                         otus.push_back(member);
1089                 }
1090
1091                 int count = 0;
1092                 for (int l = 0; l < xscores.size(); l++) {
1093
1094                         int numWithHigherRank = 0;
1095                         int numWithLowerRank = 0;
1096                         float thisrank = otus[l].score;
1097
1098                         for (int u = l+1; u < xscores.size(); u++) {
1099                                 if (otus[u].score > thisrank) { numWithHigherRank++; }
1100                                 else if (otus[u].score < thisrank) { numWithLowerRank++; }
1101                                 count++;
1102                         }
1103
1104                         numCoor += numWithHigherRank;
1105                         numDisCoor += numWithLowerRank;
1106                 }
1107
1108                 double p = (numCoor - numDisCoor) / (float) count;
1109
1110                 sig = calcKendallSig(x.size(), p);
1111
1112                 return p;
1113         }
1114         catch(exception& e) {
1115                 m->errorOut(e, "LinearAlgebra", "calcKendall");
1116                 exit(1);
1117         }
1118 }
1119 double LinearAlgebra::ran0(int& idum)
1120 {
1121     const int IA=16807,IM=2147483647,IQ=127773;
1122     const int IR=2836,MASK=123459876;
1123     const double AM=1.0/double(IM);
1124     int k;
1125     double ans;
1126
1127     idum ^= MASK;
1128     k=idum/IQ;
1129     idum=IA*(idum-k*IQ)-IR*k;
1130     if (idum < 0) idum += IM;
1131     ans=AM*idum;
1132     idum ^= MASK;
1133     return ans;
1134 }
1135
1136 double LinearAlgebra::ran1(int &idum)
1137 {
1138         const int IA=16807,IM=2147483647,IQ=127773,IR=2836,NTAB=32;
1139         const int NDIV=(1+(IM-1)/NTAB);
1140         const double EPS=3.0e-16,AM=1.0/IM,RNMX=(1.0-EPS);
1141         static int iy=0;
1142         static vector<int> iv(NTAB);
1143         int j,k;
1144         double temp;
1145
1146         if (idum <= 0 || !iy) {
1147                 if (-idum < 1) idum=1;
1148                 else idum = -idum;
1149                 for (j=NTAB+7;j>=0;j--) {
1150                         k=idum/IQ;
1151                         idum=IA*(idum-k*IQ)-IR*k;
1152                         if (idum < 0) idum += IM;
1153                         if (j < NTAB) iv[j] = idum;
1154                 }
1155                 iy=iv[0];
1156         }
1157         k=idum/IQ;
1158         idum=IA*(idum-k*IQ)-IR*k;
1159         if (idum < 0) idum += IM;
1160         j=iy/NDIV;
1161         iy=iv[j];
1162         iv[j] = idum;
1163         if ((temp=AM*iy) > RNMX) return RNMX;
1164         else return temp;
1165 }
1166
1167 double LinearAlgebra::ran2(int &idum)
1168 {
1169         const int IM1=2147483563,IM2=2147483399;
1170         const int IA1=40014,IA2=40692,IQ1=53668,IQ2=52774;
1171         const int IR1=12211,IR2=3791,NTAB=32,IMM1=IM1-1;
1172         const int NDIV=1+IMM1/NTAB;
1173         const double EPS=3.0e-16,RNMX=1.0-EPS,AM=1.0/double(IM1);
1174         static int idum2=123456789,iy=0;
1175         static vector<int> iv(NTAB);
1176         int j,k;
1177         double temp;
1178
1179         if (idum <= 0) {
1180                 idum=(idum==0 ? 1 : -idum);
1181                 idum2=idum;
1182                 for (j=NTAB+7;j>=0;j--) {
1183                         k=idum/IQ1;
1184                         idum=IA1*(idum-k*IQ1)-k*IR1;
1185                         if (idum < 0) idum += IM1;
1186                         if (j < NTAB) iv[j] = idum;
1187                 }
1188                 iy=iv[0];
1189         }
1190         k=idum/IQ1;
1191         idum=IA1*(idum-k*IQ1)-k*IR1;
1192         if (idum < 0) idum += IM1;
1193         k=idum2/IQ2;
1194         idum2=IA2*(idum2-k*IQ2)-k*IR2;
1195         if (idum2 < 0) idum2 += IM2;
1196         j=iy/NDIV;
1197         iy=iv[j]-idum2;
1198         iv[j] = idum;
1199         if (iy < 1) iy += IMM1;
1200         if ((temp=AM*iy) > RNMX) return RNMX;
1201         else return temp;
1202 }
1203
1204 double LinearAlgebra::ran3(int &idum)
1205 {
1206         static int inext,inextp;
1207         static int iff=0;
1208         const int MBIG=1000000000,MSEED=161803398,MZ=0;
1209         const double FAC=(1.0/MBIG);
1210         static vector<int> ma(56);
1211         int i,ii,k,mj,mk;
1212
1213         if (idum < 0 || iff == 0) {
1214                 iff=1;
1215                 mj=labs(MSEED-labs(idum));
1216                 mj %= MBIG;
1217                 ma[55]=mj;
1218                 mk=1;
1219                 for (i=1;i<=54;i++) {
1220                         ii=(21*i) % 55;
1221                         ma[ii]=mk;
1222                         mk=mj-mk;
1223                         if (mk < int(MZ)) mk += MBIG;
1224                         mj=ma[ii];
1225                 }
1226                 for (k=0;k<4;k++)
1227                         for (i=1;i<=55;i++) {
1228                                 ma[i] -= ma[1+(i+30) % 55];
1229                                 if (ma[i] < int(MZ)) ma[i] += MBIG;
1230                         }
1231                 inext=0;
1232                 inextp=31;
1233                 idum=1;
1234         }
1235         if (++inext == 56) inext=1;
1236         if (++inextp == 56) inextp=1;
1237         mj=ma[inext]-ma[inextp];
1238         if (mj < int(MZ)) mj += MBIG;
1239         ma[inext]=mj;
1240         return mj*FAC;
1241 }
1242
1243 double LinearAlgebra::ran4(int &idum)
1244 {
1245 #if defined(vax) || defined(_vax_) || defined(__vax__) || defined(VAX)
1246         static const unsigned long jflone = 0x00004080;
1247         static const unsigned long jflmsk = 0xffff007f;
1248 #else
1249         static const unsigned long jflone = 0x3f800000;
1250         static const unsigned long jflmsk = 0x007fffff;
1251 #endif
1252         unsigned long irword,itemp,lword;
1253         static int idums = 0;
1254
1255         if (idum < 0) {
1256                 idums = -idum;
1257                 idum=1;
1258         }
1259         irword=idum;
1260         lword=idums;
1261         psdes(lword,irword);
1262         itemp=jflone | (jflmsk & irword);
1263         ++idum;
1264         return (*(float *)&itemp)-1.0;
1265 }
1266
1267 void LinearAlgebra::psdes(unsigned long &lword, unsigned long &irword)
1268 {
1269         const int NITER=4;
1270         static const unsigned long c1[NITER]={
1271                 0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L};
1272         static const unsigned long c2[NITER]={
1273                 0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L};
1274         unsigned long i,ia,ib,iswap,itmph=0,itmpl=0;
1275
1276         for (i=0;i<NITER;i++) {
1277                 ia=(iswap=irword) ^ c1[i];
1278                 itmpl = ia & 0xffff;
1279                 itmph = ia >> 16;
1280                 ib=itmpl*itmpl+ ~(itmph*itmph);
1281                 irword=lword ^ (((ia = (ib >> 16) |
1282                           ((ib & 0xffff) << 16)) ^ c2[i])+itmpl*itmph);
1283                 lword=iswap;
1284         }
1285 }
1286 /*********************************************************************************************************************************/
1287 double LinearAlgebra::calcKendallSig(double n, double r){
1288     try {
1289
1290         double sig = 0.0;
1291         double svar=(4.0*n+10.0)/(9.0*n*(n-1.0));
1292         double z= r/sqrt(svar);
1293         sig=erfcc(fabs(z)/1.4142136);
1294
1295                 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1296
1297         return sig;
1298     }
1299         catch(exception& e) {
1300                 m->errorOut(e, "LinearAlgebra", "calcKendallSig");
1301                 exit(1);
1302         }
1303 }
1304 /*********************************************************************************************************************************/
1305 double LinearAlgebra::calcWilcoxon(vector<double>& x, vector<double>& y, double& sig){
1306         try {
1307                 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1308
1309                 double W = 0.0;
1310         sig = 0.0;
1311
1312         vector<double> signPairs;
1313         vector<spearmanRank> absV;
1314         for (int i = 0; i < x.size(); i++) {
1315             if (m->control_pressed) { return W; }
1316             double temp = y[i]-x[i];
1317             double signV = sign(temp);
1318             if (signV != 0) { // exclude zeros
1319                 spearmanRank member(toString(i), abs(temp));
1320                 absV.push_back(member);
1321                 signPairs.push_back(signV);
1322             }
1323         }
1324
1325         //rank absV
1326         //sort xscores
1327                 sort(absV.begin(), absV.end(), compareSpearman);
1328
1329                 //convert scores to ranks of x
1330                 vector<spearmanRank*> ties;
1331                 int rankTotal = 0;
1332                 for (int j = 0; j < absV.size(); j++) {
1333             if (m->control_pressed) { return W; }
1334                         rankTotal += (j+1);
1335                         ties.push_back(&(absV[j]));
1336
1337                         if (j != absV.size()-1) { // you are not the last so you can look ahead
1338                                 if (absV[j].score != absV[j+1].score) { // you are done with ties, rank them and continue
1339                                         for (int k = 0; k < ties.size(); k++) {
1340                                                 float thisrank = rankTotal / (float) ties.size();
1341                                                 (*ties[k]).score = thisrank;
1342                                         }
1343                                         ties.clear();
1344                                         rankTotal = 0;
1345                                 }
1346                         }else { // you are the last one
1347                                 for (int k = 0; k < ties.size(); k++) {
1348                                         float thisrank = rankTotal / (float) ties.size();
1349                                         (*ties[k]).score = thisrank;
1350                                 }
1351                         }
1352                 }
1353
1354         //sum ranks times sign
1355         for (int i = 0; i < absV.size(); i++) {
1356             if (m->control_pressed) { return W; }
1357             W += (absV[i].score*signPairs[i]);
1358         }
1359         W = abs(W);
1360
1361         //find zScore
1362         cout << "still need to find sig!!" << endl;
1363
1364                 return W;
1365         }
1366         catch(exception& e) {
1367                 m->errorOut(e, "LinearAlgebra", "calcWilcoxon");
1368                 exit(1);
1369         }
1370 }
1371 /*********************************************************************************************************************************/
1372 double LinearAlgebra::calcSpearman(vector<double>& x, vector<double>& y, double& sig){
1373         try {
1374                 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1375
1376                 //format data
1377         double sf = 0.0; //f^3 - f where f is the number of ties in x;
1378         double sg = 0.0; //f^3 - f where f is the number of ties in y;
1379                 map<float, int> tableX;
1380                 map<float, int>::iterator itTable;
1381                 vector<spearmanRank> xscores;
1382
1383                 for (int i = 0; i < x.size(); i++) {
1384                         spearmanRank member(toString(i), x[i]);
1385                         xscores.push_back(member);
1386
1387                         //count number of repeats
1388                         itTable = tableX.find(x[i]);
1389                         if (itTable == tableX.end()) {
1390                                 tableX[x[i]] = 1;
1391                         }else {
1392                                 tableX[x[i]]++;
1393                         }
1394                 }
1395
1396
1397                 //calc LX
1398                 double Lx = 0.0;
1399                 for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) {
1400                         double tx = (double) itTable->second;
1401                         Lx += ((pow(tx, 3.0) - tx) / 12.0);
1402                 }
1403
1404
1405                 //sort x
1406                 sort(xscores.begin(), xscores.end(), compareSpearman);
1407
1408                 //convert scores to ranks of x
1409                 //convert to ranks
1410                 map<string, float> rankx;
1411                 vector<spearmanRank> xties;
1412                 int rankTotal = 0;
1413                 for (int j = 0; j < xscores.size(); j++) {
1414                         rankTotal += (j+1);
1415                         xties.push_back(xscores[j]);
1416
1417                         if (j != xscores.size()-1) { // you are not the last so you can look ahead
1418                                 if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
1419                                         for (int k = 0; k < xties.size(); k++) {
1420                                                 float thisrank = rankTotal / (float) xties.size();
1421                                                 rankx[xties[k].name] = thisrank;
1422                                         }
1423                     int t = xties.size();
1424                     sf += (t*t*t-t);
1425                                         xties.clear();
1426                                         rankTotal = 0;
1427                                 }
1428                         }else { // you are the last one
1429                                 for (int k = 0; k < xties.size(); k++) {
1430                                         float thisrank = rankTotal / (float) xties.size();
1431                                         rankx[xties[k].name] = thisrank;
1432                                 }
1433                         }
1434                 }
1435
1436                 //format x
1437                 vector<spearmanRank> yscores;
1438                 map<float, int> tableY;
1439                 for (int j = 0; j < y.size(); j++) {
1440                         spearmanRank member(toString(j), y[j]);
1441                         yscores.push_back(member);
1442
1443                         itTable = tableY.find(member.score);
1444                         if (itTable == tableY.end()) {
1445                                 tableY[member.score] = 1;
1446                         }else {
1447                                 tableY[member.score]++;
1448                         }
1449
1450                 }
1451
1452                 //calc Ly
1453                 double Ly = 0.0;
1454                 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
1455                         double ty = (double) itTable->second;
1456                         Ly += ((pow(ty, 3.0) - ty) / 12.0);
1457                 }
1458
1459                 sort(yscores.begin(), yscores.end(), compareSpearman);
1460
1461                 //convert to ranks
1462                 map<string, float> rank;
1463                 vector<spearmanRank> yties;
1464                 rankTotal = 0;
1465                 for (int j = 0; j < yscores.size(); j++) {
1466                         rankTotal += (j+1);
1467                         yties.push_back(yscores[j]);
1468
1469                         if (j != yscores.size()-1) { // you are not the last so you can look ahead
1470                                 if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
1471                                         for (int k = 0; k < yties.size(); k++) {
1472                                                 float thisrank = rankTotal / (float) yties.size();
1473                                                 rank[yties[k].name] = thisrank;
1474                                         }
1475                     int t = yties.size();
1476                     sg += (t*t*t-t);
1477                                         yties.clear();
1478                                         rankTotal = 0;
1479                                 }
1480                         }else { // you are the last one
1481                                 for (int k = 0; k < yties.size(); k++) {
1482                                         float thisrank = rankTotal / (float) yties.size();
1483                                         rank[yties[k].name] = thisrank;
1484                                 }
1485                         }
1486                 }
1487
1488                 double di = 0.0;
1489                 for (int k = 0; k < x.size(); k++) {
1490
1491                         float xi = rankx[toString(k)];
1492                         float yi = rank[toString(k)];
1493
1494                         di += ((xi - yi) * (xi - yi));
1495                 }
1496
1497                 double p = 0.0;
1498
1499                 double n = (double) x.size();
1500                 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx;
1501                 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
1502
1503                 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
1504
1505                 //Numerical Recipes 646
1506         sig = calcSpearmanSig(n, sf, sg, di);
1507
1508                 return p;
1509         }
1510         catch(exception& e) {
1511                 m->errorOut(e, "LinearAlgebra", "calcSpearman");
1512                 exit(1);
1513         }
1514 }
1515 /*********************************************************************************************************************************/
1516 double LinearAlgebra::calcSpearmanSig(double n, double sf, double sg, double d){
1517     try {
1518
1519         double sig = 0.0;
1520         double probrs = 0.0;
1521         double en=n;
1522         double en3n=en*en*en-en;
1523         double aved=en3n/6.0-(sf+sg)/12.0;
1524         double fac=(1.0-sf/en3n)*(1.0-sg/en3n);
1525         double vard=((en-1.0)*en*en*SQR(en+1.0)/36.0)*fac;
1526         double zd=(d-aved)/sqrt(vard);
1527         double probd=erfcc(fabs(zd)/1.4142136);
1528         double rs=(1.0-(6.0/en3n)*(d+(sf+sg)/12.0))/sqrt(fac);
1529         fac=(rs+1.0)*(1.0-rs);
1530         if (fac > 0.0) {
1531             double t=rs*sqrt((en-2.0)/fac);
1532             double df=en-2.0;
1533             probrs=betai(0.5*df,0.5,df/(df+t*t));
1534         }else {
1535             probrs = 0.0;
1536         }
1537
1538         //smaller of probd and probrs is sig
1539         sig = probrs;
1540         if (probd < probrs) { sig = probd; }
1541
1542                 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1543
1544         return sig;
1545     }
1546         catch(exception& e) {
1547                 m->errorOut(e, "LinearAlgebra", "calcSpearmanSig");
1548                 exit(1);
1549         }
1550 }
1551 /*********************************************************************************************************************************/
1552 double LinearAlgebra::calcPearson(vector<double>& x, vector<double>& y, double& sig){
1553         try {
1554                 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1555
1556                 //find average X
1557                 float averageX = 0.0;
1558                 for (int i = 0; i < x.size(); i++) { averageX += x[i];  }
1559                 averageX = averageX / (float) x.size();
1560
1561                 //find average Y
1562                 float sumY = 0.0;
1563                 for (int j = 0; j < y.size(); j++) { sumY += y[j]; }
1564                 float Ybar = sumY / (float) y.size();
1565
1566                 double r = 0.0;
1567                 double numerator = 0.0;
1568                 double denomTerm1 = 0.0;
1569                 double denomTerm2 = 0.0;
1570
1571                 for (int j = 0; j < x.size(); j++) {
1572                         float Yi = y[j];
1573                         float Xi = x[j];
1574
1575                         numerator += ((Xi - averageX) * (Yi - Ybar));
1576                         denomTerm1 += ((Xi - averageX) * (Xi - averageX));
1577                         denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
1578                 }
1579
1580                 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
1581
1582                 r = numerator / denom;
1583
1584                 //Numerical Recipes pg.644
1585         sig = calcPearsonSig(x.size(), r);
1586
1587                 return r;
1588         }
1589         catch(exception& e) {
1590                 m->errorOut(e, "LinearAlgebra", "calcPearson");
1591                 exit(1);
1592         }
1593 }
1594 /*********************************************************************************************************************************/
1595 double LinearAlgebra::calcPearsonSig(double n, double r){
1596     try {
1597
1598         double sig = 0.0;
1599         const double TINY = 1.0e-20;
1600         double z = 0.5*log((1.0+r+TINY)/(1.0-r+TINY)); //Fisher's z transformation
1601
1602         //code below was giving an error in betacf with sop files
1603         //int df = n-2;
1604         //double t = r*sqrt(df/((1.0-r+TINY)*(1.0+r+TINY)));
1605         //sig = betai(0.5+df, 0.5, df/(df+t*t));
1606
1607         //Numerical Recipes says code below gives approximately the same result
1608         sig = erfcc(fabs(z*sqrt(n-1.0))/1.4142136);
1609                 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1610
1611         return sig;
1612     }
1613         catch(exception& e) {
1614                 m->errorOut(e, "LinearAlgebra", "calcPearsonSig");
1615                 exit(1);
1616         }
1617 }
1618 /*********************************************************************************************************************************/
1619
1620 vector<vector<double> > LinearAlgebra::getObservedEuclideanDistance(vector<vector<double> >& relAbundData){
1621     try {
1622
1623         int numSamples = relAbundData.size();
1624         int numOTUs = relAbundData[0].size();
1625
1626         vector<vector<double> > dMatrix(numSamples);
1627         for(int i=0;i<numSamples;i++){
1628             dMatrix[i].resize(numSamples);
1629         }
1630
1631         for(int i=0;i<numSamples;i++){
1632             for(int j=0;j<numSamples;j++){
1633
1634                 if (m->control_pressed) { return dMatrix; }
1635
1636                 double d = 0;
1637                 for(int k=0;k<numOTUs;k++){
1638                     d += pow((relAbundData[i][k] - relAbundData[j][k]), 2.0000);
1639                 }
1640                 dMatrix[i][j] = pow(d, 0.50000);
1641                 dMatrix[j][i] = dMatrix[i][j];
1642
1643             }
1644         }
1645         return dMatrix;
1646         }
1647         catch(exception& e) {
1648                 m->errorOut(e, "LinearAlgebra", "getObservedEuclideanDistance");
1649                 exit(1);
1650         }
1651 }
1652
1653 /*********************************************************************************************************************************/
1654 vector<double> LinearAlgebra::solveEquations(vector<vector<double> > A, vector<double> b){
1655     try {
1656         int length = (int)b.size();
1657         vector<double> x(length, 0);
1658         vector<int> index(length);
1659         for(int i=0;i<length;i++){  index[i] = i;   }
1660         double d;
1661
1662         ludcmp(A, index, d);  if (m->control_pressed) { return b; }
1663         lubksb(A, index, b);
1664
1665         return b;
1666     }
1667         catch(exception& e) {
1668                 m->errorOut(e, "LinearAlgebra", "solveEquations");
1669                 exit(1);
1670         }
1671 }
1672 /*********************************************************************************************************************************/
1673 vector<float> LinearAlgebra::solveEquations(vector<vector<float> > A, vector<float> b){
1674     try {
1675         int length = (int)b.size();
1676         vector<double> x(length, 0);
1677         vector<int> index(length);
1678         for(int i=0;i<length;i++){  index[i] = i;   }
1679         float d;
1680
1681         ludcmp(A, index, d);  if (m->control_pressed) { return b; }
1682         lubksb(A, index, b);
1683
1684         return b;
1685     }
1686         catch(exception& e) {
1687                 m->errorOut(e, "LinearAlgebra", "solveEquations");
1688                 exit(1);
1689         }
1690 }
1691
1692 /*********************************************************************************************************************************/
1693
1694 void LinearAlgebra::ludcmp(vector<vector<double> >& A, vector<int>& index, double& d){
1695     try {
1696         double tiny = 1e-20;
1697
1698         int n = (int)A.size();
1699         vector<double> vv(n, 0.0);
1700         double temp;
1701         int imax;
1702
1703         d = 1.0;
1704
1705         for(int i=0;i<n;i++){
1706             double big = 0.0;
1707             for(int j=0;j<n;j++){   if((temp=fabs(A[i][j])) > big ) big=temp;  }
1708             if(big==0.0){   m->mothurOut("Singular matrix in routine ludcmp\n");    }
1709             vv[i] = 1.0/big;
1710         }
1711
1712         for(int j=0;j<n;j++){
1713             if (m->control_pressed) { break; }
1714             for(int i=0;i<j;i++){
1715                 double sum = A[i][j];
1716                 for(int k=0;k<i;k++){   sum -= A[i][k] * A[k][j];   }
1717                 A[i][j] = sum;
1718             }
1719
1720             double big = 0.0;
1721             for(int i=j;i<n;i++){
1722                 double sum = A[i][j];
1723                 for(int k=0;k<j;k++){   sum -= A[i][k] * A[k][j];   }
1724                 A[i][j] = sum;
1725                 double dum;
1726                 if((dum = vv[i] * fabs(sum)) >= big){
1727                     big = dum;
1728                     imax = i;
1729                 }
1730             }
1731             if(j != imax){
1732                 for(int k=0;k<n;k++){
1733                     double dum = A[imax][k];
1734                     A[imax][k] = A[j][k];
1735                     A[j][k] = dum;
1736                 }
1737                 d = -d;
1738                 vv[imax] = vv[j];
1739             }
1740             index[j] = imax;
1741
1742             if(A[j][j] == 0.0){ A[j][j] = tiny; }
1743
1744             if(j != n-1){
1745                 double dum = 1.0/A[j][j];
1746                 for(int i=j+1;i<n;i++){ A[i][j] *= dum; }
1747             }
1748         }
1749     }
1750         catch(exception& e) {
1751                 m->errorOut(e, "LinearAlgebra", "ludcmp");
1752                 exit(1);
1753         }
1754 }
1755
1756 /*********************************************************************************************************************************/
1757
1758 void LinearAlgebra::lubksb(vector<vector<double> >& A, vector<int>& index, vector<double>& b){
1759     try {
1760         double total;
1761         int n = (int)A.size();
1762         int ii = 0;
1763
1764         for(int i=0;i<n;i++){
1765             if (m->control_pressed) { break; }
1766             int ip = index[i];
1767             total = b[ip];
1768             b[ip] = b[i];
1769
1770             if (ii != 0) {
1771                 for(int j=ii-1;j<i;j++){
1772                     total -= A[i][j] * b[j];
1773                 }
1774             }
1775             else if(total != 0){  ii = i+1;   }
1776             b[i] = total;
1777         }
1778         for(int i=n-1;i>=0;i--){
1779             total = b[i];
1780             for(int j=i+1;j<n;j++){ total -= A[i][j] * b[j];  }
1781             b[i] = total / A[i][i];
1782         }
1783     }
1784         catch(exception& e) {
1785                 m->errorOut(e, "LinearAlgebra", "lubksb");
1786                 exit(1);
1787         }
1788 }
1789 /*********************************************************************************************************************************/
1790
1791 void LinearAlgebra::ludcmp(vector<vector<float> >& A, vector<int>& index, float& d){
1792     try {
1793         double tiny = 1e-20;
1794
1795         int n = (int)A.size();
1796         vector<float> vv(n, 0.0);
1797         double temp;
1798         int imax;
1799
1800         d = 1.0;
1801
1802         for(int i=0;i<n;i++){
1803             float big = 0.0;
1804             for(int j=0;j<n;j++){   if((temp=fabs(A[i][j])) > big ) big=temp;  }
1805             if(big==0.0){   m->mothurOut("Singular matrix in routine ludcmp\n");    }
1806             vv[i] = 1.0/big;
1807         }
1808
1809         for(int j=0;j<n;j++){
1810             if (m->control_pressed) { break; }
1811             for(int i=0;i<j;i++){
1812                 float sum = A[i][j];
1813                 for(int k=0;k<i;k++){   sum -= A[i][k] * A[k][j];   }
1814                 A[i][j] = sum;
1815             }
1816
1817             float big = 0.0;
1818             for(int i=j;i<n;i++){
1819                 float sum = A[i][j];
1820                 for(int k=0;k<j;k++){   sum -= A[i][k] * A[k][j];   }
1821                 A[i][j] = sum;
1822                 float dum;
1823                 if((dum = vv[i] * fabs(sum)) >= big){
1824                     big = dum;
1825                     imax = i;
1826                 }
1827             }
1828             if(j != imax){
1829                 for(int k=0;k<n;k++){
1830                     float dum = A[imax][k];
1831                     A[imax][k] = A[j][k];
1832                     A[j][k] = dum;
1833                 }
1834                 d = -d;
1835                 vv[imax] = vv[j];
1836             }
1837             index[j] = imax;
1838
1839             if(A[j][j] == 0.0){ A[j][j] = tiny; }
1840
1841             if(j != n-1){
1842                 float dum = 1.0/A[j][j];
1843                 for(int i=j+1;i<n;i++){ A[i][j] *= dum; }
1844             }
1845         }
1846     }
1847         catch(exception& e) {
1848                 m->errorOut(e, "LinearAlgebra", "ludcmp");
1849                 exit(1);
1850         }
1851 }
1852
1853 /*********************************************************************************************************************************/
1854
1855 void LinearAlgebra::lubksb(vector<vector<float> >& A, vector<int>& index, vector<float>& b){
1856     try {
1857         float total;
1858         int n = (int)A.size();
1859         int ii = 0;
1860
1861         for(int i=0;i<n;i++){
1862             if (m->control_pressed) { break; }
1863             int ip = index[i];
1864             total = b[ip];
1865             b[ip] = b[i];
1866
1867             if (ii != 0) {
1868                 for(int j=ii-1;j<i;j++){
1869                     total -= A[i][j] * b[j];
1870                 }
1871             }
1872             else if(total != 0){  ii = i+1;   }
1873             b[i] = total;
1874         }
1875         for(int i=n-1;i>=0;i--){
1876             total = b[i];
1877             for(int j=i+1;j<n;j++){ total -= A[i][j] * b[j];  }
1878             b[i] = total / A[i][i];
1879         }
1880     }
1881         catch(exception& e) {
1882                 m->errorOut(e, "LinearAlgebra", "lubksb");
1883                 exit(1);
1884         }
1885 }
1886
1887
1888 /*********************************************************************************************************************************/
1889
1890 vector<vector<double> > LinearAlgebra::getInverse(vector<vector<double> > matrix){
1891     try {
1892         int n = (int)matrix.size();
1893
1894         vector<vector<double> > inverse(n);
1895         for(int i=0;i<n;i++){   inverse[i].assign(n, 0.0000);   }
1896
1897         vector<double> column(n, 0.0000);
1898         vector<int> index(n, 0);
1899         double dummy;
1900
1901         ludcmp(matrix, index, dummy);
1902
1903         for(int j=0;j<n;j++){
1904             if (m->control_pressed) { break; }
1905
1906             column.assign(n, 0);
1907
1908             column[j] = 1.0000;
1909
1910             lubksb(matrix, index, column);
1911
1912             for(int i=0;i<n;i++){   inverse[i][j] = column[i];  }
1913         }
1914
1915         return inverse;
1916     }
1917         catch(exception& e) {
1918                 m->errorOut(e, "LinearAlgebra", "getInverse");
1919                 exit(1);
1920         }
1921 }
1922
1923 /*********************************************************************************************************************************/