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 #include "wilcox.h"
  12
  13 #define PI 3.1415926535897932384626433832795
  14
  15 // This class references functions used from "Numerical Recipes in C++" //
  16
  17 /*********************************************************************************************************************************/
  18 inline double SQR(const double a)
  19 {
  20     return a*a;
  21 }
  22 /*********************************************************************************************************************************/
  23
  24 inline double SIGN(const double a, const double b)
  25 {
  26     return b>=0 ? (a>=0 ? a:-a) : (a>=0 ? -a:a);
  27 }
  28 /*********************************************************************************************************************************/
  29 //NUmerical recipes pg. 245 - Returns the complementary error function erfc(x) with fractional error everywhere less than 1.2 × 10−7.
  30 double LinearAlgebra::erfcc(double x){
  31     try {
  32         double t,z,ans;
  33         z=fabs(x);
  34         t=1.0/(1.0+0.5*z);
  35
  36         ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
  37             t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
  38             t*(-0.82215223+t*0.17087277)))))))));
  39
  40         //cout << "in erfcc " << t << '\t' << ans<< endl;
  41         return (x >= 0.0 ? ans : 2.0 - ans);
  42     }
  43         catch(exception& e) {
  44                 m->errorOut(e, "LinearAlgebra", "betai");
  45                 exit(1);
  46         }
  47 }
  48 /*********************************************************************************************************************************/
  49 //Numerical Recipes pg. 232
  50 double LinearAlgebra::betai(const double a, const double b, const double x) {
  51     try {
  52         double bt;
  53         double result = 0.0;
  54
  55         if (x < 0.0 || x > 1.0) { m->mothurOut("[ERROR]: bad x in betai.\n"); m->control_pressed = true; return 0.0; }
  56
  57         if (x == 0.0 || x == 1.0)  { bt = 0.0; }
  58         else { bt = exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));  }
  59
  60         if (x < (a+1.0) / (a+b+2.0)) { result = bt*betacf(a,b,x)/a; }
  61         else { result = 1.0-bt*betacf(b,a,1.0-x)/b; }
  62
  63         return result;
  64     }
  65         catch(exception& e) {
  66                 m->errorOut(e, "LinearAlgebra", "betai");
  67                 exit(1);
  68         }
  69 }
  70 /*********************************************************************************************************************************/
  71 //Numerical Recipes pg. 219
  72 double LinearAlgebra::gammln(const double xx) {
  73     try {
  74         int j;
  75         double x,y,tmp,ser;
  76         static const double cof[6]={76.18009172947146,-86.50532032941677,24.01409824083091,
  77             -1.231739572450155,0.120858003e-2,-0.536382e-5};
  78
  79         y=x=xx;
  80         tmp=x+5.5;
  81         tmp -= (x+0.5)*log(tmp);
  82         ser=1.0;
  83         for (j=0;j<6;j++) {
  84             ser += cof[j]/++y;
  85         }
  86         return -tmp+log(2.5066282746310005*ser/x);
  87     }
  88         catch(exception& e) {
  89                 m->errorOut(e, "LinearAlgebra", "gammln");
  90                 exit(1);
  91         }
  92 }
  93 /*********************************************************************************************************************************/
  94 //Numerical Recipes pg. 223
  95 double LinearAlgebra::gammp(const double a, const double x) {
  96     try {
  97         double gamser,gammcf,gln;
  98
  99         if (x < 0.0 || a <= 0.0) { m->mothurOut("[ERROR]: Invalid arguments in routine GAMMP\n"); m->control_pressed = true; return 0.0;}
 100         if (x < (a+1.0)) {
 101             gser(gamser,a,x,gln);
 102             return gamser;
 103         } else {
 104             gcf(gammcf,a,x,gln);
 105             return 1.0-gammcf;
 106         }
 107
 108         return 0;
 109     }
 110         catch(exception& e) {
 111                 m->errorOut(e, "LinearAlgebra", "gammp");
 112                 exit(1);
 113         }
 114 }
 115 /*********************************************************************************************************************************/
 116 //Numerical Recipes pg. 223
 117 /*double LinearAlgebra::gammq(const double a, const double x) {
 118     try {
 119         double gamser,gammcf,gln;
 120
 121         if (x < 0.0 || a <= 0.0) { m->mothurOut("[ERROR]: Invalid arguments in routine GAMMQ\n"); m->control_pressed = true; return 0.0; }
 122         if (x < (a+1.0)) {
 123             gser(gamser,a,x,gln);
 124             return 1.0-gamser;
 125         } else {
 126             gcf(gammcf,a,x,gln);
 127             return gammcf;
 128         }
 129
 130         return 0;
 131     }
 132         catch(exception& e) {
 133                 m->errorOut(e, "LinearAlgebra", "gammq");
 134                 exit(1);
 135         }
 136 }
 137 *********************************************************************************************************************************/
 138 //Numerical Recipes pg. 224
 139 double LinearAlgebra::gcf(double& gammcf, const double a, const double x, double& gln){
 140     try {
 141         const int ITMAX=100;
 142         const double EPS=numeric_limits<double>::epsilon();
 143         const double FPMIN=numeric_limits<double>::min()/EPS;
 144         int i;
 145         double an,b,c,d,del,h;
 146
 147         gln=gammln(a);
 148         b=x+1.0-a;
 149         c=1.0/FPMIN;
 150         d=1.0/b;
 151         h=d;
 152         for (i=1;i<=ITMAX;i++) {
 153             an = -i*(i-a);
 154             b += 2.0;
 155             d=an*d+b;
 156             if (fabs(d) < FPMIN) { d=FPMIN; }
 157             c=b+an/c;
 158             if (fabs(c) < FPMIN) { c=FPMIN; }
 159             d=1.0/d;
 160             del=d*c;
 161             h *= del;
 162             if (fabs(del-1.0) <= EPS) break;
 163         }
 164         if (i > ITMAX)  { m->mothurOut("[ERROR]: " + toString(a) + " too large, ITMAX=100 too small in gcf\n"); m->control_pressed = true; }
 165         gammcf=exp(-x+a*log(x)-gln)*h;
 166
 167         return 0.0;
 168     }
 169         catch(exception& e) {
 170                 m->errorOut(e, "LinearAlgebra", "gcf");
 171                 exit(1);
 172         }
 173
 174 }
 175 /*********************************************************************************************************************************/
 176 //Numerical Recipes pg. 223
 177 double LinearAlgebra::gser(double& gamser, const double a, const double x, double& gln) {
 178     try {
 179         int n;
 180         double sum,del,ap;
 181         const double EPS = numeric_limits<double>::epsilon();
 182
 183         gln=gammln(a);
 184         if (x <= 0.0) {
 185             if (x < 0.0) {  m->mothurOut("[ERROR]: x less than 0 in routine GSER\n"); m->control_pressed = true;  }
 186             gamser=0.0; return 0.0;
 187         } else {
 188             ap=a;
 189             del=sum=1.0/a;
 190             for (n=0;n<100;n++) {
 191                 ++ap;
 192                 del *= x/ap;
 193                 sum += del;
 194                 if (fabs(del) < fabs(sum)*EPS) {
 195                     gamser=sum*exp(-x+a*log(x)-gln);
 196                     return 0.0;
 197                 }
 198             }
 199
 200             m->mothurOut("[ERROR]: a too large, ITMAX too small in routine GSER\n");
 201             return 0.0;
 202         }
 203         return 0;
 204     }
 205         catch(exception& e) {
 206                 m->errorOut(e, "LinearAlgebra", "gser");
 207                 exit(1);
 208         }
 209 }
 210 /*********************************************************************************************************************************/
 211 //Numerical Recipes pg. 233
 212 double LinearAlgebra::betacf(const double a, const double b, const double x) {
 213     try {
 214         const int MAXIT = 100;
 215         const double EPS = numeric_limits<double>::epsilon();
 216         const double FPMIN = numeric_limits<double>::min() / EPS;
 217         int m1, m2;
 218         double aa, c, d, del, h, qab, qam, qap;
 219
 220         qab=a+b;
 221         qap=a+1.0;
 222         qam=a-1.0;
 223         c=1.0;
 224         d=1.0-qab*x/qap;
 225         if (fabs(d) < FPMIN) d=FPMIN;
 226         d=1.0/d;
 227         h=d;
 228         for (m1=1;m1<=MAXIT;m1++) {
 229             m2=2*m1;
 230             aa=m1*(b-m1)*x/((qam+m2)*(a+m2));
 231             d=1.0+aa*d;
 232             if (fabs(d) < FPMIN) d=FPMIN;
 233             c=1.0+aa/c;
 234             if (fabs(c) < FPMIN) c=FPMIN;
 235             d=1.0/d;
 236             h *= d*c;
 237             aa = -(a+m1)*(qab+m1)*x/((a+m2)*(qap+m2));
 238             d=1.0+aa*d;
 239             if (fabs(d) < FPMIN) d=FPMIN;
 240             c=1.0+aa/c;
 241             if (fabs(c) < FPMIN) c=FPMIN;
 242             d=1.0/d;
 243             del=d*c;
 244             h *= del;
 245             if (fabs(del-1.0) < EPS) break;
 246         }
 247
 248         if (m1 > MAXIT) { m->mothurOut("[ERROR]: a or b too big or MAXIT too small in betacf."); m->mothurOutEndLine(); m->control_pressed = true; }
 249         return h;
 250
 251     }
 252         catch(exception& e) {
 253                 m->errorOut(e, "LinearAlgebra", "betacf");
 254                 exit(1);
 255         }
 256 }
 257 /*********************************************************************************************************************************/
 258 //[3][4] * [4][5] - columns in first must match rows in second, returns matrix[3][5]
 259 vector<vector<double> > LinearAlgebra::matrix_mult(vector<vector<double> > first, vector<vector<double> > second){
 260         try {
 261                 vector<vector<double> > product;
 262
 263                 int first_rows = first.size();
 264                 int first_cols = first[0].size();
 265                 int second_cols = second[0].size();
 266
 267                 product.resize(first_rows);
 268                 for(int i=0;i<first_rows;i++){
 269                         product[i].resize(second_cols);
 270                 }
 271
 272                 for(int i=0;i<first_rows;i++){
 273                         for(int j=0;j<second_cols;j++){
 274
 275                                 if (m->control_pressed) { return product; }
 276
 277                                 product[i][j] = 0.0;
 278                                 for(int k=0;k<first_cols;k++){
 279                                         product[i][j] += first[i][k] * second[k][j];
 280                                 }
 281                         }
 282                 }
 283
 284                 return product;
 285         }
 286         catch(exception& e) {
 287                 m->errorOut(e, "LinearAlgebra", "matrix_mult");
 288                 exit(1);
 289         }
 290
 291 }
 292 /*********************************************************************************************************************************/
 293
 294 vector<vector<double> > LinearAlgebra::transpose(vector<vector<double> >matrix){
 295         try {
 296                 vector<vector<double> > trans; trans.resize(matrix[0].size());
 297         for (int i = 0; i < trans.size(); i++) {
 298             for (int j = 0; j < matrix.size(); j++) { trans[i].push_back(matrix[j][i]); }
 299         }
 300
 301                 return trans;
 302         }
 303         catch(exception& e) {
 304                 m->errorOut(e, "LinearAlgebra", "transpose");
 305                 exit(1);
 306         }
 307
 308 }
 309 /*********************************************************************************************************************************/
 310
 311 void LinearAlgebra::recenter(double offset, vector<vector<double> > D, vector<vector<double> >& G){
 312         try {
 313                 int rank = D.size();
 314
 315                 vector<vector<double> > A(rank);
 316                 vector<vector<double> > C(rank);
 317                 for(int i=0;i<rank;i++){
 318                         A[i].resize(rank);
 319                         C[i].resize(rank);
 320                 }
 321
 322                 double scale = -1.0000 / (double) rank;
 323
 324                 for(int i=0;i<rank;i++){
 325                         A[i][i] = 0.0000;
 326                         C[i][i] = 1.0000 + scale;
 327                         for(int j=i+1;j<rank;j++){
 328                                 A[i][j] = A[j][i] = -0.5 * D[i][j] * D[i][j] + offset;
 329                                 C[i][j] = C[j][i] = scale;
 330                         }
 331                 }
 332
 333                 A = matrix_mult(C,A);
 334                 G = matrix_mult(A,C);
 335         }
 336         catch(exception& e) {
 337                 m->errorOut(e, "LinearAlgebra", "recenter");
 338                 exit(1);
 339         }
 340
 341 }
 342 /*********************************************************************************************************************************/
 343
 344 //  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
 345
 346 int LinearAlgebra::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
 347         try {
 348                 double scale, hh, h, g, f;
 349
 350                 int n = a.size();
 351
 352                 d.resize(n);
 353                 e.resize(n);
 354
 355                 for(int i=n-1;i>0;i--){
 356                         int l=i-1;
 357                         h = scale = 0.0000;
 358                         if(l>0){
 359                                 for(int k=0;k<l+1;k++){
 360                                         scale += fabs(a[i][k]);
 361                                 }
 362                                 if(scale == 0.0){
 363                                         e[i] = a[i][l];
 364                                 }
 365                                 else{
 366                                         for(int k=0;k<l+1;k++){
 367                                                 a[i][k] /= scale;
 368                                                 h += a[i][k] * a[i][k];
 369                                         }
 370                                         f = a[i][l];
 371                                         g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
 372                                         e[i] = scale * g;
 373                                         h -= f * g;
 374                                         a[i][l] = f - g;
 375                                         f = 0.0;
 376                                         for(int j=0;j<l+1;j++){
 377                                                 a[j][i] = a[i][j] / h;
 378                                                 g = 0.0;
 379                                                 for(int k=0;k<j+1;k++){
 380                                                         g += a[j][k] * a[i][k];
 381                                                 }
 382                                                 for(int k=j+1;k<l+1;k++){
 383                                                         g += a[k][j] * a[i][k];
 384                                                 }
 385                                                 e[j] = g / h;
 386                                                 f += e[j] * a[i][j];
 387                                         }
 388                                         hh = f / (h + h);
 389                                         for(int j=0;j<l+1;j++){
 390                                                 f = a[i][j];
 391                                                 e[j] = g = e[j] - hh * f;
 392                                                 for(int k=0;k<j+1;k++){
 393                                                         a[j][k] -= (f * e[k] + g * a[i][k]);
 394                                                 }
 395                                         }
 396                                 }
 397                         }
 398                         else{
 399                                 e[i] = a[i][l];
 400                         }
 401
 402                         d[i] = h;
 403                 }
 404
 405                 d[0] = 0.0000;
 406                 e[0] = 0.0000;
 407
 408                 for(int i=0;i<n;i++){
 409                         int l = i;
 410                         if(d[i] != 0.0){
 411                                 for(int j=0;j<l;j++){
 412                                         g = 0.0000;
 413                                         for(int k=0;k<l;k++){
 414                                                 g += a[i][k] * a[k][j];
 415                                         }
 416                                         for(int k=0;k<l;k++){
 417                                                 a[k][j] -= g * a[k][i];
 418                                         }
 419                                 }
 420                         }
 421                         d[i] = a[i][i];
 422                         a[i][i] = 1.0000;
 423                         for(int j=0;j<l;j++){
 424                                 a[j][i] = a[i][j] = 0.0;
 425                         }
 426                 }
 427
 428                 return 0;
 429         }
 430         catch(exception& e) {
 431                 m->errorOut(e, "LinearAlgebra", "tred2");
 432                 exit(1);
 433         }
 434
 435 }
 436 /*********************************************************************************************************************************/
 437
 438 double LinearAlgebra::pythag(double a, double b)        {       return(pow(a*a+b*b,0.5));       }
 439
 440 /*********************************************************************************************************************************/
 441
 442 //  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
 443
 444 int LinearAlgebra::qtli(vector<double>& d, vector<double>& e, vector<vector<double> >& z) {
 445         try {
 446                 int myM, i, iter;
 447                 double s, r, p, g, f, dd, c, b;
 448
 449                 int n = d.size();
 450                 for(int i=1;i<=n;i++){
 451                         e[i-1] = e[i];
 452                 }
 453                 e[n-1] = 0.0000;
 454
 455                 for(int l=0;l<n;l++){
 456                         iter = 0;
 457                         do {
 458                                 for(myM=l;myM<n-1;myM++){
 459                                         dd = fabs(d[myM]) + fabs(d[myM+1]);
 460                                         if(fabs(e[myM])+dd == dd) break;
 461                                 }
 462                                 if(myM != l){
 463                                         if(iter++ == 3000) cerr << "Too many iterations in tqli\n";
 464                                         g = (d[l+1]-d[l]) / (2.0 * e[l]);
 465                                         r = pythag(g, 1.0);
 466                                         g = d[myM] - d[l] + e[l] / (g + SIGN(r,g));
 467                                         s = c = 1.0;
 468                                         p = 0.0000;
 469                                         for(i=myM-1;i>=l;i--){
 470                                                 f = s * e[i];
 471                                                 b = c * e[i];
 472                                                 e[i+1] = (r=pythag(f,g));
 473                                                 if(r==0.0){
 474                                                         d[i+1] -= p;
 475                                                         e[myM] = 0.0000;
 476                                                         break;
 477                                                 }
 478                                                 s = f / r;
 479                                                 c = g / r;
 480                                                 g = d[i+1] - p;
 481                                                 r = (d[i] - g) * s + 2.0 * c * b;
 482                                                 d[i+1] = g + ( p = s * r);
 483                                                 g = c * r - b;
 484                                                 for(int k=0;k<n;k++){
 485                                                         f = z[k][i+1];
 486                                                         z[k][i+1] = s * z[k][i] + c * f;
 487                                                         z[k][i] = c * z[k][i] - s * f;
 488                                                 }
 489                                         }
 490                                         if(r == 0.00 && i >= l) continue;
 491                                         d[l] -= p;
 492                                         e[l] = g;
 493                                         e[myM] = 0.0;
 494                                 }
 495                         } while (myM != l);
 496                 }
 497
 498                 int k;
 499                 for(int i=0;i<n;i++){
 500                         p=d[k=i];
 501                         for(int j=i;j<n;j++){
 502                                 if(d[j] >= p){
 503                                         p=d[k=j];
 504                                 }
 505                         }
 506                         if(k!=i){
 507                                 d[k]=d[i];
 508                                 d[i]=p;
 509                                 for(int j=0;j<n;j++){
 510                                         p=z[j][i];
 511                                         z[j][i] = z[j][k];
 512                                         z[j][k] = p;
 513                                 }
 514                         }
 515                 }
 516
 517                 return 0;
 518         }
 519         catch(exception& e) {
 520                 m->errorOut(e, "LinearAlgebra", "qtli");
 521                 exit(1);
 522         }
 523 }
 524 /*********************************************************************************************************************************/
 525 //groups by dimension
 526 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes, int dimensions){
 527         try {
 528                 //make square matrix
 529                 vector< vector<double> > dists; dists.resize(axes.size());
 530                 for (int i = 0; i < dists.size(); i++) {  dists[i].resize(axes.size(), 0.0); }
 531
 532                 if (dimensions == 1) { //one dimension calc = abs(x-y)
 533
 534                         for (int i = 0; i < dists.size(); i++) {
 535
 536                                 if (m->control_pressed) { return dists; }
 537
 538                                 for (int j = 0; j < i; j++) {
 539                                         dists[i][j] = abs(axes[i][0] - axes[j][0]);
 540                                         dists[j][i] = dists[i][j];
 541                                 }
 542                         }
 543
 544                 }else if (dimensions > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
 545
 546                         for (int i = 0; i < dists.size(); i++) {
 547
 548                                 if (m->control_pressed) { return dists; }
 549
 550                                 for (int j = 0; j < i; j++) {
 551                                         double sum = 0.0;
 552                                         for (int k = 0; k < dimensions; k++) {
 553                                                 sum += ((axes[i][k] - axes[j][k]) * (axes[i][k] - axes[j][k]));
 554                                         }
 555
 556                                         dists[i][j] = sqrt(sum);
 557                                         dists[j][i] = dists[i][j];
 558                                 }
 559                         }
 560
 561                 }
 562
 563                 return dists;
 564         }
 565         catch(exception& e) {
 566                 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
 567                 exit(1);
 568         }
 569 }
 570 /*********************************************************************************************************************************/
 571 //returns groups by dimensions from dimensions by groups
 572 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes){
 573         try {
 574                 //make square matrix
 575                 vector< vector<double> > dists; dists.resize(axes[0].size());
 576                 for (int i = 0; i < dists.size(); i++) {  dists[i].resize(axes[0].size(), 0.0); }
 577
 578                 if (axes.size() == 1) { //one dimension calc = abs(x-y)
 579
 580                         for (int i = 0; i < dists.size(); i++) {
 581
 582                                 if (m->control_pressed) { return dists; }
 583
 584                                 for (int j = 0; j < i; j++) {
 585                                         dists[i][j] = abs(axes[0][i] - axes[0][j]);
 586                                         dists[j][i] = dists[i][j];
 587                                 }
 588                         }
 589
 590                 }else if (axes.size() > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
 591
 592                         for (int i = 0; i < dists[0].size(); i++) {
 593
 594                                 if (m->control_pressed) { return dists; }
 595
 596                                 for (int j = 0; j < i; j++) {
 597                                         double sum = 0.0;
 598                                         for (int k = 0; k < axes.size(); k++) {
 599                                                 sum += ((axes[k][i] - axes[k][j]) * (axes[k][i] - axes[k][j]));
 600                                         }
 601
 602                                         dists[i][j] = sqrt(sum);
 603                                         dists[j][i] = dists[i][j];
 604                                 }
 605                         }
 606
 607                 }
 608
 609                 return dists;
 610         }
 611         catch(exception& e) {
 612                 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
 613                 exit(1);
 614         }
 615 }
 616 /*********************************************************************************************************************************/
 617 //assumes both matrices are square and the same size
 618 double LinearAlgebra::calcPearson(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
 619         try {
 620
 621                 //find average for - X
 622                 int count = 0;
 623                 float averageEuclid = 0.0;
 624                 for (int i = 0; i < euclidDists.size(); i++) {
 625                         for (int j = 0; j < i; j++) {
 626                                 averageEuclid += euclidDists[i][j];
 627                                 count++;
 628                         }
 629                 }
 630                 averageEuclid = averageEuclid / (float) count;
 631
 632                 //find average for - Y
 633                 count = 0;
 634                 float averageUser = 0.0;
 635                 for (int i = 0; i < userDists.size(); i++) {
 636                         for (int j = 0; j < i; j++) {
 637                                 averageUser += userDists[i][j];
 638                                 count++;
 639                         }
 640                 }
 641                 averageUser = averageUser / (float) count;
 642
 643                 double numerator = 0.0;
 644                 double denomTerm1 = 0.0;
 645                 double denomTerm2 = 0.0;
 646
 647                 for (int i = 0; i < euclidDists.size(); i++) {
 648
 649                         for (int k = 0; k < i; k++) { //just lt dists
 650
 651                                 float Yi = userDists[i][k];
 652                                 float Xi = euclidDists[i][k];
 653
 654                                 numerator += ((Xi - averageEuclid) * (Yi - averageUser));
 655                                 denomTerm1 += ((Xi - averageEuclid) * (Xi - averageEuclid));
 656                                 denomTerm2 += ((Yi - averageUser) * (Yi - averageUser));
 657                         }
 658                 }
 659
 660                 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
 661                 double r = numerator / denom;
 662
 663                 //divide by zero error
 664                 if (isnan(r) || isinf(r)) { r = 0.0; }
 665
 666                 return r;
 667
 668         }
 669         catch(exception& e) {
 670                 m->errorOut(e, "LinearAlgebra", "calcPearson");
 671                 exit(1);
 672         }
 673 }
 674 /*********************************************************************************************************************************/
 675 //assumes both matrices are square and the same size
 676 double LinearAlgebra::calcSpearman(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
 677         try {
 678                 double r;
 679
 680                 //format data
 681                 map<float, int> tableX;
 682                 map<float, int>::iterator itTable;
 683                 vector<spearmanRank> scores;
 684
 685                 for (int i = 0; i < euclidDists.size(); i++) {
 686                         for (int j = 0; j < i; j++) {
 687                                 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
 688                                 scores.push_back(member);
 689
 690                                 //count number of repeats
 691                                 itTable = tableX.find(euclidDists[i][j]);
 692                                 if (itTable == tableX.end()) {
 693                                         tableX[euclidDists[i][j]] = 1;
 694                                 }else {
 695                                         tableX[euclidDists[i][j]]++;
 696                                 }
 697                         }
 698                 }
 699
 700                 //sort scores
 701                 sort(scores.begin(), scores.end(), compareSpearman);
 702
 703                 //calc LX
 704                 double Lx = 0.0;
 705                 for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) {
 706                         double tx = (double) itTable->second;
 707                         Lx += ((pow(tx, 3.0) - tx) / 12.0);
 708                 }
 709
 710                 //find ranks of xi
 711                 map<string, float> rankEuclid;
 712                 vector<spearmanRank> ties;
 713                 int rankTotal = 0;
 714                 for (int j = 0; j < scores.size(); j++) {
 715                         rankTotal += (j+1);
 716                         ties.push_back(scores[j]);
 717
 718                         if (j != (scores.size()-1)) { // you are not the last so you can look ahead
 719                                 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
 720
 721                                         for (int k = 0; k < ties.size(); k++) {
 722                                                 float thisrank = rankTotal / (float) ties.size();
 723                                                 rankEuclid[ties[k].name] = thisrank;
 724                                         }
 725                                         ties.clear();
 726                                         rankTotal = 0;
 727                                 }
 728                         }else { // you are the last one
 729
 730                                 for (int k = 0; k < ties.size(); k++) {
 731                                         float thisrank = rankTotal / (float) ties.size();
 732                                         rankEuclid[ties[k].name] = thisrank;
 733                                 }
 734                         }
 735                 }
 736
 737
 738                 //format data
 739                 map<float, int> tableY;
 740                 scores.clear();
 741
 742                 for (int i = 0; i < userDists.size(); i++) {
 743                         for (int j = 0; j < i; j++) {
 744                                 spearmanRank member(toString(scores.size()), userDists[i][j]);
 745                                 scores.push_back(member);
 746
 747                                 //count number of repeats
 748                                 itTable = tableY.find(userDists[i][j]);
 749                                 if (itTable == tableY.end()) {
 750                                         tableY[userDists[i][j]] = 1;
 751                                 }else {
 752                                         tableY[userDists[i][j]]++;
 753                                 }
 754                         }
 755                 }
 756
 757                 //sort scores
 758                 sort(scores.begin(), scores.end(), compareSpearman);
 759
 760                 //calc LX
 761                 double Ly = 0.0;
 762                 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
 763                         double ty = (double) itTable->second;
 764                         Ly += ((pow(ty, 3.0) - ty) / 12.0);
 765                 }
 766
 767                 //find ranks of yi
 768                 map<string, float> rankUser;
 769                 ties.clear();
 770                 rankTotal = 0;
 771                 for (int j = 0; j < scores.size(); j++) {
 772                         rankTotal += (j+1);
 773                         ties.push_back(scores[j]);
 774
 775                         if (j != (scores.size()-1)) { // you are not the last so you can look ahead
 776                                 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
 777
 778                                         for (int k = 0; k < ties.size(); k++) {
 779                                                 float thisrank = rankTotal / (float) ties.size();
 780                                                 rankUser[ties[k].name] = thisrank;
 781                                         }
 782                                         ties.clear();
 783                                         rankTotal = 0;
 784                                 }
 785                         }else { // you are the last one
 786
 787                                 for (int k = 0; k < ties.size(); k++) {
 788                                         float thisrank = rankTotal / (float) ties.size();
 789                                         rankUser[ties[k].name] = thisrank;
 790                                 }
 791                         }
 792                 }
 793
 794
 795                 double di = 0.0;
 796                 int count = 0;
 797                 for (int i = 0; i < userDists.size(); i++) {
 798                         for (int j = 0; j < i; j++) {
 799
 800                                 float xi = rankEuclid[toString(count)];
 801                                 float yi = rankUser[toString(count)];
 802
 803                                 di += ((xi - yi) * (xi - yi));
 804
 805                                 count++;
 806                         }
 807                 }
 808
 809                 double n = (double) count;
 810
 811                 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx;
 812                 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
 813
 814                 r = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
 815
 816                 //divide by zero error
 817                 if (isnan(r) || isinf(r)) { r = 0.0; }
 818
 819                 return r;
 820
 821         }
 822         catch(exception& e) {
 823                 m->errorOut(e, "LinearAlgebra", "calcSpearman");
 824                 exit(1);
 825         }
 826 }
 827
 828 /*********************************************************************************************************************************/
 829 //assumes both matrices are square and the same size
 830 double LinearAlgebra::calcKendall(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
 831         try {
 832                 double r;
 833
 834                 //format data
 835                 vector<spearmanRank> scores;
 836                 for (int i = 0; i < euclidDists.size(); i++) {
 837                         for (int j = 0; j < i; j++) {
 838                                 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
 839                                 scores.push_back(member);
 840                         }
 841                 }
 842
 843                 //sort scores
 844                 sort(scores.begin(), scores.end(), compareSpearman);
 845
 846                 //find ranks of xi
 847                 map<string, float> rankEuclid;
 848                 vector<spearmanRank> ties;
 849                 int rankTotal = 0;
 850                 for (int j = 0; j < scores.size(); j++) {
 851                         rankTotal += (j+1);
 852                         ties.push_back(scores[j]);
 853
 854                         if (j != (scores.size()-1)) { // you are not the last so you can look ahead
 855                                 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
 856
 857                                         for (int k = 0; k < ties.size(); k++) {
 858                                                 float thisrank = rankTotal / (float) ties.size();
 859                                                 rankEuclid[ties[k].name] = thisrank;
 860                                         }
 861                                         ties.clear();
 862                                         rankTotal = 0;
 863                                 }
 864                         }else { // you are the last one
 865
 866                                 for (int k = 0; k < ties.size(); k++) {
 867                                         float thisrank = rankTotal / (float) ties.size();
 868                                         rankEuclid[ties[k].name] = thisrank;
 869                                 }
 870                         }
 871                 }
 872
 873                 vector<spearmanRank> scoresUser;
 874                 for (int i = 0; i < userDists.size(); i++) {
 875                         for (int j = 0; j < i; j++) {
 876                                 spearmanRank member(toString(scoresUser.size()), userDists[i][j]);
 877                                 scoresUser.push_back(member);
 878                         }
 879                 }
 880
 881                 //sort scores
 882                 sort(scoresUser.begin(), scoresUser.end(), compareSpearman);
 883
 884                 //find ranks of yi
 885                 map<string, float> rankUser;
 886                 ties.clear();
 887                 rankTotal = 0;
 888                 for (int j = 0; j < scoresUser.size(); j++) {
 889                         rankTotal += (j+1);
 890                         ties.push_back(scoresUser[j]);
 891
 892                         if (j != (scoresUser.size()-1)) { // you are not the last so you can look ahead
 893                                 if (scoresUser[j].score != scoresUser[j+1].score) { // you are done with ties, rank them and continue
 894
 895                                         for (int k = 0; k < ties.size(); k++) {
 896                                                 float thisrank = rankTotal / (float) ties.size();
 897                                                 rankUser[ties[k].name] = thisrank;
 898                                         }
 899                                         ties.clear();
 900                                         rankTotal = 0;
 901                                 }
 902                         }else { // you are the last one
 903
 904                                 for (int k = 0; k < ties.size(); k++) {
 905                                         float thisrank = rankTotal / (float) ties.size();
 906                                         rankUser[ties[k].name] = thisrank;
 907                                 }
 908                         }
 909                 }
 910
 911                 int numCoor = 0;
 912                 int numDisCoor = 0;
 913
 914                 //order user ranks
 915                 vector<spearmanRank> user;
 916                 for (int l = 0; l < scores.size(); l++) {
 917                         spearmanRank member(scores[l].name, rankUser[scores[l].name]);
 918                         user.push_back(member);
 919                 }
 920
 921                 int count = 0;
 922                 for (int l = 0; l < scores.size(); l++) {
 923
 924                         int numWithHigherRank = 0;
 925                         int numWithLowerRank = 0;
 926                         float thisrank = user[l].score;
 927
 928                         for (int u = l+1; u < scores.size(); u++) {
 929                                 if (user[u].score > thisrank) { numWithHigherRank++; }
 930                                 else if (user[u].score < thisrank) { numWithLowerRank++; }
 931                                 count++;
 932                         }
 933
 934                         numCoor += numWithHigherRank;
 935                         numDisCoor += numWithLowerRank;
 936                 }
 937
 938                 r = (numCoor - numDisCoor) / (float) count;
 939
 940                 //divide by zero error
 941                 if (isnan(r) || isinf(r)) { r = 0.0; }
 942
 943                 return r;
 944
 945         }
 946         catch(exception& e) {
 947                 m->errorOut(e, "LinearAlgebra", "calcKendall");
 948                 exit(1);
 949         }
 950 }
 951 /*********************************************************************************************************************************/
 952 double LinearAlgebra::calcKendall(vector<double>& x, vector<double>& y, double& sig){
 953         try {
 954                 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
 955
 956                 //format data
 957                 vector<spearmanRank> xscores;
 958                 for (int i = 0; i < x.size(); i++) {
 959                         spearmanRank member(toString(i), x[i]);
 960                         xscores.push_back(member);
 961                 }
 962
 963                 //sort xscores
 964                 sort(xscores.begin(), xscores.end(), compareSpearman);
 965
 966                 //convert scores to ranks of x
 967                 vector<spearmanRank*> ties;
 968                 int rankTotal = 0;
 969                 for (int j = 0; j < xscores.size(); j++) {
 970                         rankTotal += (j+1);
 971                         ties.push_back(&(xscores[j]));
 972
 973                         if (j != xscores.size()-1) { // you are not the last so you can look ahead
 974                                 if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
 975                                         for (int k = 0; k < ties.size(); k++) {
 976                                                 float thisrank = rankTotal / (float) ties.size();
 977                                                 (*ties[k]).score = thisrank;
 978                                         }
 979                                         ties.clear();
 980                                         rankTotal = 0;
 981                                 }
 982                         }else { // you are the last one
 983                                 for (int k = 0; k < ties.size(); k++) {
 984                                         float thisrank = rankTotal / (float) ties.size();
 985                                         (*ties[k]).score = thisrank;
 986                                 }
 987                         }
 988                 }
 989
 990
 991                 //format data
 992                 vector<spearmanRank> yscores;
 993                 for (int j = 0; j < y.size(); j++) {
 994                         spearmanRank member(toString(j), y[j]);
 995                         yscores.push_back(member);
 996                 }
 997
 998                 //sort yscores
 999                 sort(yscores.begin(), yscores.end(), compareSpearman);
1000
1001                 //convert to ranks
1002                 map<string, float> rank;
1003                 vector<spearmanRank> yties;
1004                 rankTotal = 0;
1005                 for (int j = 0; j < yscores.size(); j++) {
1006                         rankTotal += (j+1);
1007                         yties.push_back(yscores[j]);
1008
1009                         if (j != yscores.size()-1) { // you are not the last so you can look ahead
1010                                 if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
1011                                         for (int k = 0; k < yties.size(); k++) {
1012                                                 float thisrank = rankTotal / (float) yties.size();
1013                                                 rank[yties[k].name] = thisrank;
1014                                         }
1015                                         yties.clear();
1016                                         rankTotal = 0;
1017                                 }
1018                         }else { // you are the last one
1019                                 for (int k = 0; k < yties.size(); k++) {
1020                                         float thisrank = rankTotal / (float) yties.size();
1021                                         rank[yties[k].name] = thisrank;
1022                                 }
1023                         }
1024                 }
1025
1026
1027                 int numCoor = 0;
1028                 int numDisCoor = 0;
1029
1030                 //associate x and y
1031                 vector<spearmanRank> otus;
1032                 for (int l = 0; l < xscores.size(); l++) {
1033                         spearmanRank member(xscores[l].name, rank[xscores[l].name]);
1034                         otus.push_back(member);
1035                 }
1036
1037                 int count = 0;
1038                 for (int l = 0; l < xscores.size(); l++) {
1039
1040                         int numWithHigherRank = 0;
1041                         int numWithLowerRank = 0;
1042                         float thisrank = otus[l].score;
1043
1044                         for (int u = l+1; u < xscores.size(); u++) {
1045                                 if (otus[u].score > thisrank) { numWithHigherRank++; }
1046                                 else if (otus[u].score < thisrank) { numWithLowerRank++; }
1047                                 count++;
1048                         }
1049
1050                         numCoor += numWithHigherRank;
1051                         numDisCoor += numWithLowerRank;
1052                 }
1053
1054                 double p = (numCoor - numDisCoor) / (float) count;
1055
1056                 sig = calcKendallSig(x.size(), p);
1057
1058                 return p;
1059         }
1060         catch(exception& e) {
1061                 m->errorOut(e, "LinearAlgebra", "calcKendall");
1062                 exit(1);
1063         }
1064 }
1065 double LinearAlgebra::ran0(int& idum)
1066 {
1067     const int IA=16807,IM=2147483647,IQ=127773;
1068     const int IR=2836,MASK=123459876;
1069     const double AM=1.0/double(IM);
1070     int k;
1071     double ans;
1072
1073     idum ^= MASK;
1074     k=idum/IQ;
1075     idum=IA*(idum-k*IQ)-IR*k;
1076     if (idum < 0) idum += IM;
1077     ans=AM*idum;
1078     idum ^= MASK;
1079     return ans;
1080 }
1081
1082 double LinearAlgebra::ran1(int &idum)
1083 {
1084         const int IA=16807,IM=2147483647,IQ=127773,IR=2836,NTAB=32;
1085         const int NDIV=(1+(IM-1)/NTAB);
1086         const double EPS=3.0e-16,AM=1.0/IM,RNMX=(1.0-EPS);
1087         static int iy=0;
1088         static vector<int> iv(NTAB);
1089         int j,k;
1090         double temp;
1091
1092         if (idum <= 0 || !iy) {
1093                 if (-idum < 1) idum=1;
1094                 else idum = -idum;
1095                 for (j=NTAB+7;j>=0;j--) {
1096                         k=idum/IQ;
1097                         idum=IA*(idum-k*IQ)-IR*k;
1098                         if (idum < 0) idum += IM;
1099                         if (j < NTAB) iv[j] = idum;
1100                 }
1101                 iy=iv[0];
1102         }
1103         k=idum/IQ;
1104         idum=IA*(idum-k*IQ)-IR*k;
1105         if (idum < 0) idum += IM;
1106         j=iy/NDIV;
1107         iy=iv[j];
1108         iv[j] = idum;
1109         if ((temp=AM*iy) > RNMX) return RNMX;
1110         else return temp;
1111 }
1112
1113 double LinearAlgebra::ran2(int &idum)
1114 {
1115         const int IM1=2147483563,IM2=2147483399;
1116         const int IA1=40014,IA2=40692,IQ1=53668,IQ2=52774;
1117         const int IR1=12211,IR2=3791,NTAB=32,IMM1=IM1-1;
1118         const int NDIV=1+IMM1/NTAB;
1119         const double EPS=3.0e-16,RNMX=1.0-EPS,AM=1.0/double(IM1);
1120         static int idum2=123456789,iy=0;
1121         static vector<int> iv(NTAB);
1122         int j,k;
1123         double temp;
1124
1125         if (idum <= 0) {
1126                 idum=(idum==0 ? 1 : -idum);
1127                 idum2=idum;
1128                 for (j=NTAB+7;j>=0;j--) {
1129                         k=idum/IQ1;
1130                         idum=IA1*(idum-k*IQ1)-k*IR1;
1131                         if (idum < 0) idum += IM1;
1132                         if (j < NTAB) iv[j] = idum;
1133                 }
1134                 iy=iv[0];
1135         }
1136         k=idum/IQ1;
1137         idum=IA1*(idum-k*IQ1)-k*IR1;
1138         if (idum < 0) idum += IM1;
1139         k=idum2/IQ2;
1140         idum2=IA2*(idum2-k*IQ2)-k*IR2;
1141         if (idum2 < 0) idum2 += IM2;
1142         j=iy/NDIV;
1143         iy=iv[j]-idum2;
1144         iv[j] = idum;
1145         if (iy < 1) iy += IMM1;
1146         if ((temp=AM*iy) > RNMX) return RNMX;
1147         else return temp;
1148 }
1149
1150 double LinearAlgebra::ran3(int &idum)
1151 {
1152         static int inext,inextp;
1153         static int iff=0;
1154         const int MBIG=1000000000,MSEED=161803398,MZ=0;
1155         const double FAC=(1.0/MBIG);
1156         static vector<int> ma(56);
1157         int i,ii,k,mj,mk;
1158
1159         if (idum < 0 || iff == 0) {
1160                 iff=1;
1161                 mj=labs(MSEED-labs(idum));
1162                 mj %= MBIG;
1163                 ma[55]=mj;
1164                 mk=1;
1165                 for (i=1;i<=54;i++) {
1166                         ii=(21*i) % 55;
1167                         ma[ii]=mk;
1168                         mk=mj-mk;
1169                         if (mk < int(MZ)) mk += MBIG;
1170                         mj=ma[ii];
1171                 }
1172                 for (k=0;k<4;k++)
1173                         for (i=1;i<=55;i++) {
1174                                 ma[i] -= ma[1+(i+30) % 55];
1175                                 if (ma[i] < int(MZ)) ma[i] += MBIG;
1176                         }
1177                 inext=0;
1178                 inextp=31;
1179                 idum=1;
1180         }
1181         if (++inext == 56) inext=1;
1182         if (++inextp == 56) inextp=1;
1183         mj=ma[inext]-ma[inextp];
1184         if (mj < int(MZ)) mj += MBIG;
1185         ma[inext]=mj;
1186         return mj*FAC;
1187 }
1188
1189 double LinearAlgebra::ran4(int &idum)
1190 {
1191 #if defined(vax) || defined(_vax_) || defined(__vax__) || defined(VAX)
1192         static const unsigned long jflone = 0x00004080;
1193         static const unsigned long jflmsk = 0xffff007f;
1194 #else
1195         static const unsigned long jflone = 0x3f800000;
1196         static const unsigned long jflmsk = 0x007fffff;
1197 #endif
1198         unsigned long irword,itemp,lword;
1199         static int idums = 0;
1200
1201         if (idum < 0) {
1202                 idums = -idum;
1203                 idum=1;
1204         }
1205         irword=idum;
1206         lword=idums;
1207         psdes(lword,irword);
1208         itemp=jflone | (jflmsk & irword);
1209         ++idum;
1210         return (*(float *)&itemp)-1.0;
1211 }
1212
1213 void LinearAlgebra::psdes(unsigned long &lword, unsigned long &irword)
1214 {
1215         const int NITER=4;
1216         static const unsigned long c1[NITER]={
1217                 0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L};
1218         static const unsigned long c2[NITER]={
1219                 0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L};
1220         unsigned long i,ia,ib,iswap,itmph=0,itmpl=0;
1221
1222         for (i=0;i<NITER;i++) {
1223                 ia=(iswap=irword) ^ c1[i];
1224                 itmpl = ia & 0xffff;
1225                 itmph = ia >> 16;
1226                 ib=itmpl*itmpl+ ~(itmph*itmph);
1227                 irword=lword ^ (((ia = (ib >> 16) |
1228                           ((ib & 0xffff) << 16)) ^ c2[i])+itmpl*itmph);
1229                 lword=iswap;
1230         }
1231 }
1232 /*********************************************************************************************************************************/
1233 double LinearAlgebra::calcKendallSig(double n, double r){
1234     try {
1235
1236         double sig = 0.0;
1237         double svar=(4.0*n+10.0)/(9.0*n*(n-1.0));
1238         double z= r/sqrt(svar);
1239         sig=erfcc(fabs(z)/1.4142136);
1240
1241                 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1242
1243         return sig;
1244     }
1245         catch(exception& e) {
1246                 m->errorOut(e, "LinearAlgebra", "calcKendallSig");
1247                 exit(1);
1248         }
1249 }
1250 /*********************************************************************************************************************************/
1251 double LinearAlgebra::calcKruskalWallis(vector<spearmanRank>& values, double& pValue){
1252         try {
1253         double H;
1254         set<string> treatments;
1255
1256         //rank values
1257         sort(values.begin(), values.end(), compareSpearman);
1258         vector<spearmanRank*> ties;
1259         int rankTotal = 0;
1260         vector<int> TIES;
1261         for (int j = 0; j < values.size(); j++) {
1262             treatments.insert(values[j].name);
1263             rankTotal += (j+1);
1264             ties.push_back(&(values[j]));
1265
1266             if (j != values.size()-1) { // you are not the last so you can look ahead
1267                 if (values[j].score != values[j+1].score) { // you are done with ties, rank them and continue
1268                     if (ties.size() > 1) { TIES.push_back(ties.size()); }
1269                     for (int k = 0; k < ties.size(); k++) {
1270                         double thisrank = rankTotal / (double) ties.size();
1271                         (*ties[k]).score = thisrank;
1272                     }
1273                     ties.clear();
1274                     rankTotal = 0;
1275                 }
1276             }else { // you are the last one
1277                 if (ties.size() > 1) { TIES.push_back(ties.size()); }
1278                 for (int k = 0; k < ties.size(); k++) {
1279                     double thisrank = rankTotal / (double) ties.size();
1280                     (*ties[k]).score = thisrank;
1281                 }
1282             }
1283         }
1284
1285
1286         // H = 12/(N*(N+1)) * (sum Ti^2/n) - 3(N+1)
1287         map<string, double> sums;
1288         map<string, double> counts;
1289         for (set<string>::iterator it = treatments.begin(); it != treatments.end(); it++) { sums[*it] = 0.0; counts[*it] = 0; }
1290
1291         for (int j = 0; j < values.size(); j++) {
1292             sums[values[j].name] += values[j].score;
1293             counts[values[j].name]+= 1.0;
1294         }
1295
1296         double middleTerm = 0.0;
1297         for (set<string>::iterator it = treatments.begin(); it != treatments.end(); it++) {
1298             middleTerm += ((sums[*it]*sums[*it])/counts[*it]);
1299         }
1300
1301         double firstTerm = 12 / (double) (values.size()*(values.size()+1));
1302         double lastTerm = 3 * (values.size()+1);
1303
1304         H = firstTerm * middleTerm - lastTerm;
1305
1306         //adjust for ties
1307         if (TIES.size() != 0) {
1308             double sum = 0.0;
1309             for (int j = 0; j < TIES.size(); j++) { sum += ((TIES[j]*TIES[j]*TIES[j])-TIES[j]); }
1310             double result = 1.0 - (sum / (double) ((values.size()*values.size()*values.size())-values.size()));
1311             H /= result;
1312         }
1313
1314         if (isnan(H) || isinf(H)) { H = 0; }
1315
1316         //Numerical Recipes pg221
1317         pValue = 1.0 - (gammp(((treatments.size()-1)/(double)2.0), H/2.0));
1318
1319         return H;
1320     }
1321         catch(exception& e) {
1322                 m->errorOut(e, "LinearAlgebra", "calcKruskalWallis");
1323                 exit(1);
1324         }
1325 }
1326 /*********************************************************************************************************************************/
1327 double LinearAlgebra::normalvariate(double mean, double standardDeviation) {
1328     try {
1329         double u1 = ((double)(rand()) + 1.0 )/( (double)(RAND_MAX) + 1.0);
1330         double u2 = ((double)(rand()) + 1.0 )/( (double)(RAND_MAX) + 1.0);
1331         //double r = sqrt( -2.0*log(u1) );
1332         //double theta = 2.0*PI*u2;
1333         //cout << cos(8.*atan(1.)*u2)*sqrt(-2.*log(u1)) << endl;
1334         return cos(8.*atan(1.)*u2)*sqrt(-2.*log(u1));
1335     }
1336         catch(exception& e) {
1337                 m->errorOut(e, "LinearAlgebra", "normalvariate");
1338                 exit(1);
1339         }
1340 }
1341 /*********************************************************************************************************************************/
1342 //thanks http://www.johndcook.com/cpp_phi.html
1343 double LinearAlgebra::pnorm(double x){
1344     try {
1345         // constants
1346         double a1 =  0.254829592;
1347         double a2 = -0.284496736;
1348         double a3 =  1.421413741;
1349         double a4 = -1.453152027;
1350         double a5 =  1.061405429;
1351         double p  =  0.3275911;
1352
1353         // Save the sign of x
1354         int sign = 1;
1355         if (x < 0)
1356             sign = -1;
1357         x = fabs(x)/sqrt(2.0);
1358
1359         // A&S formula 7.1.26
1360         double t = 1.0/(1.0 + p*x);
1361         double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);
1362
1363         return 0.5*(1.0 + sign*y);
1364
1365     }
1366         catch(exception& e) {
1367                 m->errorOut(e, "LinearAlgebra", "pnorm");
1368                 exit(1);
1369         }
1370 }
1371
1372 /*********************************************************************************************************************************/
1373 double LinearAlgebra::calcWilcoxon(vector<double>& x, vector<double>& y, double& sig){
1374         try {
1375                 double W = 0.0;
1376         sig = 0.0;
1377
1378         vector<spearmanRank> ranks;
1379         for (int i = 0; i < x.size(); i++) {
1380             if (m->control_pressed) { return W; }
1381             spearmanRank member("x", x[i]);
1382             ranks.push_back(member);
1383         }
1384
1385         for (int i = 0; i < y.size(); i++) {
1386             if (m->control_pressed) { return W; }
1387             spearmanRank member("y", y[i]);
1388             ranks.push_back(member);
1389         }
1390
1391         //sort values
1392                 sort(ranks.begin(), ranks.end(), compareSpearman);
1393
1394                 //convert scores to ranks of x
1395                 vector<spearmanRank*> ties;
1396                 int rankTotal = 0;
1397         vector<int> TIES;
1398                 for (int j = 0; j < ranks.size(); j++) {
1399             if (m->control_pressed) { return W; }
1400                         rankTotal += (j+1);
1401                         ties.push_back(&(ranks[j]));
1402
1403                         if (j != ranks.size()-1) { // you are not the last so you can look ahead
1404                                 if (ranks[j].score != ranks[j+1].score) { // you are done with ties, rank them and continue
1405                     if (ties.size() > 1) { TIES.push_back(ties.size()); }
1406                                         for (int k = 0; k < ties.size(); k++) {
1407                                                 float thisrank = rankTotal / (float) ties.size();
1408                                                 (*ties[k]).score = thisrank;
1409                                         }
1410                                         ties.clear();
1411                                         rankTotal = 0;
1412                                 }
1413                         }else { // you are the last one
1414                 if (ties.size() > 1) { TIES.push_back(ties.size()); }
1415                                 for (int k = 0; k < ties.size(); k++) {
1416                                         float thisrank = rankTotal / (float) ties.size();
1417                                         (*ties[k]).score = thisrank;
1418                                 }
1419                         }
1420                 }
1421
1422         //from R wilcox.test function
1423         //STATISTIC <- sum(r[seq_along(x)]) - n.x * (n.x + 1)/2
1424         double sumRanks = 0.0;
1425         for (int i = 0; i < ranks.size(); i++) {
1426             if (m->control_pressed) { return W; }
1427             if (ranks[i].name == "x") { sumRanks += ranks[i].score; }
1428         }
1429
1430         W = sumRanks - x.size() * ((double)(x.size() + 1)) / 2.0;
1431
1432         //exact <- (n.x < 50) && (n.y < 50)
1433         bool findExact = false;
1434         if ((x.size() < 50) && (y.size() < 50)) { findExact = true; }
1435
1436
1437         if (findExact && (TIES.size() == 0)) { //find exact and no ties
1438             //PVAL <- switch(alternative, two.sided = {
1439             //p <- if (STATISTIC > (n.x * n.y/2))
1440             PWilcox wilcox;
1441             double pval = 0.0;
1442             if (W > ((double)x.size()*y.size()/2.0)) {
1443                 //pwilcox(STATISTIC-1, n.x, n.y, lower.tail = FALSE)
1444                 pval = wilcox.pwilcox(W-1, x.size(), y.size(), false);
1445             }else {
1446                 //pwilcox(STATISTIC,n.x, n.y)
1447                 pval = wilcox.pwilcox(W, x.size(), y.size(), true);
1448             }
1449             sig = 2.0 * pval;
1450             if (1.0 < sig) { sig = 1.0; }
1451         }else {
1452             //z <- STATISTIC - n.x * n.y/2
1453             double z = W - (double)(x.size() * y.size()/2.0);
1454             //NTIES <- table(r)
1455             double sum = 0.0;
1456             for (int j = 0; j < TIES.size(); j++) { sum += ((TIES[j]*TIES[j]*TIES[j])-TIES[j]); }
1457
1458             //SIGMA <- sqrt((n.x * n.y/12) * ((n.x + n.y + 1) -
1459                                             //sum(NTIES^3 - NTIES)/((n.x + n.y) * (n.x + n.y -
1460                                                                             //1))))
1461             double sigma = 0.0;
1462             double firstTerm = (double)(x.size() * y.size()/12.0);
1463             double secondTerm = (double)(x.size() + y.size() + 1) - sum / (double)((x.size() + y.size()) * (x.size() + y.size() - 1));
1464             sigma = sqrt(firstTerm * secondTerm);
1465
1466             //CORRECTION <- switch(alternative, two.sided = sign(z) * 0.5, greater = 0.5, less = -0.5)
1467             double CORRECTION = 0.0;
1468             if (z < 0) { CORRECTION = -1.0; }
1469             else if (z > 0) { CORRECTION = 1.0; }
1470             CORRECTION *= 0.5;
1471
1472             z = (z - CORRECTION)/sigma;
1473
1474             //PVAL <- switch(alternative,  two.sided = 2 * min(pnorm(z), pnorm(z, lower.tail = FALSE)))
1475             sig = pnorm(z);
1476             if ((1.0-sig) < sig) { sig = 1.0 - sig; }
1477             sig *= 2;
1478         }
1479
1480         return W;
1481         }
1482         catch(exception& e) {
1483                 m->errorOut(e, "LinearAlgebra", "calcWilcoxon");
1484                 exit(1);
1485         }
1486 }
1487
1488 /*********************************************************************************************************************************/
1489 double LinearAlgebra::choose(double n, double k){
1490         try {
1491         n = floor(n + 0.5);
1492         k = floor(k + 0.5);
1493
1494         double lchoose = gammln(n + 1.0) - gammln(k + 1.0) - gammln(n - k + 1.0);
1495
1496         return (floor(exp(lchoose) + 0.5));
1497     }
1498         catch(exception& e) {
1499                 m->errorOut(e, "LinearAlgebra", "choose");
1500                 exit(1);
1501         }
1502 }
1503 /*********************************************************************************************************************************/
1504 double LinearAlgebra::calcSpearman(vector<double>& x, vector<double>& y, double& sig){
1505         try {
1506                 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1507
1508                 //format data
1509         double sf = 0.0; //f^3 - f where f is the number of ties in x;
1510         double sg = 0.0; //f^3 - f where f is the number of ties in y;
1511                 map<float, int> tableX;
1512                 map<float, int>::iterator itTable;
1513                 vector<spearmanRank> xscores;
1514
1515                 for (int i = 0; i < x.size(); i++) {
1516                         spearmanRank member(toString(i), x[i]);
1517                         xscores.push_back(member);
1518
1519                         //count number of repeats
1520                         itTable = tableX.find(x[i]);
1521                         if (itTable == tableX.end()) {
1522                                 tableX[x[i]] = 1;
1523                         }else {
1524                                 tableX[x[i]]++;
1525                         }
1526                 }
1527
1528
1529                 //calc LX
1530                 double Lx = 0.0;
1531                 for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) {
1532                         double tx = (double) itTable->second;
1533                         Lx += ((pow(tx, 3.0) - tx) / 12.0);
1534                 }
1535
1536
1537                 //sort x
1538                 sort(xscores.begin(), xscores.end(), compareSpearman);
1539
1540                 //convert scores to ranks of x
1541                 //convert to ranks
1542                 map<string, float> rankx;
1543                 vector<spearmanRank> xties;
1544                 int rankTotal = 0;
1545                 for (int j = 0; j < xscores.size(); j++) {
1546                         rankTotal += (j+1);
1547                         xties.push_back(xscores[j]);
1548
1549                         if (j != xscores.size()-1) { // you are not the last so you can look ahead
1550                                 if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
1551                                         for (int k = 0; k < xties.size(); k++) {
1552                                                 float thisrank = rankTotal / (float) xties.size();
1553                                                 rankx[xties[k].name] = thisrank;
1554                                         }
1555                     int t = xties.size();
1556                     sf += (t*t*t-t);
1557                                         xties.clear();
1558                                         rankTotal = 0;
1559                                 }
1560                         }else { // you are the last one
1561                                 for (int k = 0; k < xties.size(); k++) {
1562                                         float thisrank = rankTotal / (float) xties.size();
1563                                         rankx[xties[k].name] = thisrank;
1564                                 }
1565                         }
1566                 }
1567
1568                 //format x
1569                 vector<spearmanRank> yscores;
1570                 map<float, int> tableY;
1571                 for (int j = 0; j < y.size(); j++) {
1572                         spearmanRank member(toString(j), y[j]);
1573                         yscores.push_back(member);
1574
1575                         itTable = tableY.find(member.score);
1576                         if (itTable == tableY.end()) {
1577                                 tableY[member.score] = 1;
1578                         }else {
1579                                 tableY[member.score]++;
1580                         }
1581
1582                 }
1583
1584                 //calc Ly
1585                 double Ly = 0.0;
1586                 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
1587                         double ty = (double) itTable->second;
1588                         Ly += ((pow(ty, 3.0) - ty) / 12.0);
1589                 }
1590
1591                 sort(yscores.begin(), yscores.end(), compareSpearman);
1592
1593                 //convert to ranks
1594                 map<string, float> rank;
1595                 vector<spearmanRank> yties;
1596                 rankTotal = 0;
1597                 for (int j = 0; j < yscores.size(); j++) {
1598                         rankTotal += (j+1);
1599                         yties.push_back(yscores[j]);
1600
1601                         if (j != yscores.size()-1) { // you are not the last so you can look ahead
1602                                 if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
1603                                         for (int k = 0; k < yties.size(); k++) {
1604                                                 float thisrank = rankTotal / (float) yties.size();
1605                                                 rank[yties[k].name] = thisrank;
1606                                         }
1607                     int t = yties.size();
1608                     sg += (t*t*t-t);
1609                                         yties.clear();
1610                                         rankTotal = 0;
1611                                 }
1612                         }else { // you are the last one
1613                                 for (int k = 0; k < yties.size(); k++) {
1614                                         float thisrank = rankTotal / (float) yties.size();
1615                                         rank[yties[k].name] = thisrank;
1616                                 }
1617                         }
1618                 }
1619
1620                 double di = 0.0;
1621                 for (int k = 0; k < x.size(); k++) {
1622
1623                         float xi = rankx[toString(k)];
1624                         float yi = rank[toString(k)];
1625
1626                         di += ((xi - yi) * (xi - yi));
1627                 }
1628
1629                 double p = 0.0;
1630
1631                 double n = (double) x.size();
1632                 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx;
1633                 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
1634
1635                 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
1636
1637                 //Numerical Recipes 646
1638         sig = calcSpearmanSig(n, sf, sg, di);
1639
1640                 return p;
1641         }
1642         catch(exception& e) {
1643                 m->errorOut(e, "LinearAlgebra", "calcSpearman");
1644                 exit(1);
1645         }
1646 }
1647 /*********************************************************************************************************************************/
1648 double LinearAlgebra::calcSpearmanSig(double n, double sf, double sg, double d){
1649     try {
1650
1651         double sig = 0.0;
1652         double probrs = 0.0;
1653         double en=n;
1654         double en3n=en*en*en-en;
1655         double aved=en3n/6.0-(sf+sg)/12.0;
1656         double fac=(1.0-sf/en3n)*(1.0-sg/en3n);
1657         double vard=((en-1.0)*en*en*SQR(en+1.0)/36.0)*fac;
1658         double zd=(d-aved)/sqrt(vard);
1659         double probd=erfcc(fabs(zd)/1.4142136);
1660         double rs=(1.0-(6.0/en3n)*(d+(sf+sg)/12.0))/sqrt(fac);
1661         fac=(rs+1.0)*(1.0-rs);
1662         if (fac > 0.0) {
1663             double t=rs*sqrt((en-2.0)/fac);
1664             double df=en-2.0;
1665             probrs=betai(0.5*df,0.5,df/(df+t*t));
1666         }else {
1667             probrs = 0.0;
1668         }
1669
1670         //smaller of probd and probrs is sig
1671         sig = probrs;
1672         if (probd < probrs) { sig = probd; }
1673
1674                 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1675
1676         return sig;
1677     }
1678         catch(exception& e) {
1679                 m->errorOut(e, "LinearAlgebra", "calcSpearmanSig");
1680                 exit(1);
1681         }
1682 }
1683 /*********************************************************************************************************************************/
1684 double LinearAlgebra::calcPearson(vector<double>& x, vector<double>& y, double& sig){
1685         try {
1686                 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1687
1688                 //find average X
1689                 float averageX = 0.0;
1690                 for (int i = 0; i < x.size(); i++) { averageX += x[i];  }
1691                 averageX = averageX / (float) x.size();
1692
1693                 //find average Y
1694                 float sumY = 0.0;
1695                 for (int j = 0; j < y.size(); j++) { sumY += y[j]; }
1696                 float Ybar = sumY / (float) y.size();
1697
1698                 double r = 0.0;
1699                 double numerator = 0.0;
1700                 double denomTerm1 = 0.0;
1701                 double denomTerm2 = 0.0;
1702
1703                 for (int j = 0; j < x.size(); j++) {
1704                         float Yi = y[j];
1705                         float Xi = x[j];
1706
1707                         numerator += ((Xi - averageX) * (Yi - Ybar));
1708                         denomTerm1 += ((Xi - averageX) * (Xi - averageX));
1709                         denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
1710                 }
1711
1712                 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
1713
1714                 r = numerator / denom;
1715
1716                 //Numerical Recipes pg.644
1717         sig = calcPearsonSig(x.size(), r);
1718
1719                 return r;
1720         }
1721         catch(exception& e) {
1722                 m->errorOut(e, "LinearAlgebra", "calcPearson");
1723                 exit(1);
1724         }
1725 }
1726 /*********************************************************************************************************************************/
1727 double LinearAlgebra::calcPearsonSig(double n, double r){
1728     try {
1729
1730         double sig = 0.0;
1731         const double TINY = 1.0e-20;
1732         double z = 0.5*log((1.0+r+TINY)/(1.0-r+TINY)); //Fisher's z transformation
1733
1734         //code below was giving an error in betacf with sop files
1735         //int df = n-2;
1736         //double t = r*sqrt(df/((1.0-r+TINY)*(1.0+r+TINY)));
1737         //sig = betai(0.5+df, 0.5, df/(df+t*t));
1738
1739         //Numerical Recipes says code below gives approximately the same result
1740         sig = erfcc(fabs(z*sqrt(n-1.0))/1.4142136);
1741                 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1742
1743         return sig;
1744     }
1745         catch(exception& e) {
1746                 m->errorOut(e, "LinearAlgebra", "calcPearsonSig");
1747                 exit(1);
1748         }
1749 }
1750 /*********************************************************************************************************************************/
1751
1752 vector<vector<double> > LinearAlgebra::getObservedEuclideanDistance(vector<vector<double> >& relAbundData){
1753     try {
1754
1755         int numSamples = relAbundData.size();
1756         int numOTUs = relAbundData[0].size();
1757
1758         vector<vector<double> > dMatrix(numSamples);
1759         for(int i=0;i<numSamples;i++){
1760             dMatrix[i].resize(numSamples);
1761         }
1762
1763         for(int i=0;i<numSamples;i++){
1764             for(int j=0;j<numSamples;j++){
1765
1766                 if (m->control_pressed) { return dMatrix; }
1767
1768                 double d = 0;
1769                 for(int k=0;k<numOTUs;k++){
1770                     d += pow((relAbundData[i][k] - relAbundData[j][k]), 2.0000);
1771                 }
1772                 dMatrix[i][j] = pow(d, 0.50000);
1773                 dMatrix[j][i] = dMatrix[i][j];
1774
1775             }
1776         }
1777         return dMatrix;
1778         }
1779         catch(exception& e) {
1780                 m->errorOut(e, "LinearAlgebra", "getObservedEuclideanDistance");
1781                 exit(1);
1782         }
1783 }
1784
1785 /*********************************************************************************************************************************/
1786 vector<double> LinearAlgebra::solveEquations(vector<vector<double> > A, vector<double> b){
1787     try {
1788         int length = (int)b.size();
1789         vector<double> x(length, 0);
1790         vector<int> index(length);
1791         for(int i=0;i<length;i++){  index[i] = i;   }
1792         double d;
1793
1794         ludcmp(A, index, d);  if (m->control_pressed) { return b; }
1795         lubksb(A, index, b);
1796
1797         return b;
1798     }
1799         catch(exception& e) {
1800                 m->errorOut(e, "LinearAlgebra", "solveEquations");
1801                 exit(1);
1802         }
1803 }
1804 /*********************************************************************************************************************************/
1805 vector<float> LinearAlgebra::solveEquations(vector<vector<float> > A, vector<float> b){
1806     try {
1807         int length = (int)b.size();
1808         vector<double> x(length, 0);
1809         vector<int> index(length);
1810         for(int i=0;i<length;i++){  index[i] = i;   }
1811         float d;
1812
1813         ludcmp(A, index, d);  if (m->control_pressed) { return b; }
1814         lubksb(A, index, b);
1815
1816         return b;
1817     }
1818         catch(exception& e) {
1819                 m->errorOut(e, "LinearAlgebra", "solveEquations");
1820                 exit(1);
1821         }
1822 }
1823
1824 /*********************************************************************************************************************************/
1825
1826 void LinearAlgebra::ludcmp(vector<vector<double> >& A, vector<int>& index, double& d){
1827     try {
1828         double tiny = 1e-20;
1829
1830         int n = (int)A.size();
1831         vector<double> vv(n, 0.0);
1832         double temp;
1833         int imax;
1834
1835         d = 1.0;
1836
1837         for(int i=0;i<n;i++){
1838             double big = 0.0;
1839             for(int j=0;j<n;j++){   if((temp=fabs(A[i][j])) > big ) big=temp;  }
1840             if(big==0.0){   m->mothurOut("Singular matrix in routine ludcmp\n");    }
1841             vv[i] = 1.0/big;
1842         }
1843
1844         for(int j=0;j<n;j++){
1845             if (m->control_pressed) { break; }
1846             for(int i=0;i<j;i++){
1847                 double sum = A[i][j];
1848                 for(int k=0;k<i;k++){   sum -= A[i][k] * A[k][j];   }
1849                 A[i][j] = sum;
1850             }
1851
1852             double big = 0.0;
1853             for(int i=j;i<n;i++){
1854                 double sum = A[i][j];
1855                 for(int k=0;k<j;k++){   sum -= A[i][k] * A[k][j];   }
1856                 A[i][j] = sum;
1857                 double dum;
1858                 if((dum = vv[i] * fabs(sum)) >= big){
1859                     big = dum;
1860                     imax = i;
1861                 }
1862             }
1863             if(j != imax){
1864                 for(int k=0;k<n;k++){
1865                     double dum = A[imax][k];
1866                     A[imax][k] = A[j][k];
1867                     A[j][k] = dum;
1868                 }
1869                 d = -d;
1870                 vv[imax] = vv[j];
1871             }
1872             index[j] = imax;
1873
1874             if(A[j][j] == 0.0){ A[j][j] = tiny; }
1875
1876             if(j != n-1){
1877                 double dum = 1.0/A[j][j];
1878                 for(int i=j+1;i<n;i++){ A[i][j] *= dum; }
1879             }
1880         }
1881     }
1882         catch(exception& e) {
1883                 m->errorOut(e, "LinearAlgebra", "ludcmp");
1884                 exit(1);
1885         }
1886 }
1887
1888 /*********************************************************************************************************************************/
1889
1890 void LinearAlgebra::lubksb(vector<vector<double> >& A, vector<int>& index, vector<double>& b){
1891     try {
1892         double total;
1893         int n = (int)A.size();
1894         int ii = 0;
1895
1896         for(int i=0;i<n;i++){
1897             if (m->control_pressed) { break; }
1898             int ip = index[i];
1899             total = b[ip];
1900             b[ip] = b[i];
1901
1902             if (ii != 0) {
1903                 for(int j=ii-1;j<i;j++){
1904                     total -= A[i][j] * b[j];
1905                 }
1906             }
1907             else if(total != 0){  ii = i+1;   }
1908             b[i] = total;
1909         }
1910         for(int i=n-1;i>=0;i--){
1911             total = b[i];
1912             for(int j=i+1;j<n;j++){ total -= A[i][j] * b[j];  }
1913             b[i] = total / A[i][i];
1914         }
1915     }
1916         catch(exception& e) {
1917                 m->errorOut(e, "LinearAlgebra", "lubksb");
1918                 exit(1);
1919         }
1920 }
1921 /*********************************************************************************************************************************/
1922
1923 void LinearAlgebra::ludcmp(vector<vector<float> >& A, vector<int>& index, float& d){
1924     try {
1925         double tiny = 1e-20;
1926
1927         int n = (int)A.size();
1928         vector<float> vv(n, 0.0);
1929         double temp;
1930         int imax;
1931
1932         d = 1.0;
1933
1934         for(int i=0;i<n;i++){
1935             float big = 0.0;
1936             for(int j=0;j<n;j++){   if((temp=fabs(A[i][j])) > big ) big=temp;  }
1937             if(big==0.0){   m->mothurOut("Singular matrix in routine ludcmp\n");    }
1938             vv[i] = 1.0/big;
1939         }
1940
1941         for(int j=0;j<n;j++){
1942             if (m->control_pressed) { break; }
1943             for(int i=0;i<j;i++){
1944                 float sum = A[i][j];
1945                 for(int k=0;k<i;k++){   sum -= A[i][k] * A[k][j];   }
1946                 A[i][j] = sum;
1947             }
1948
1949             float big = 0.0;
1950             for(int i=j;i<n;i++){
1951                 float sum = A[i][j];
1952                 for(int k=0;k<j;k++){   sum -= A[i][k] * A[k][j];   }
1953                 A[i][j] = sum;
1954                 float dum;
1955                 if((dum = vv[i] * fabs(sum)) >= big){
1956                     big = dum;
1957                     imax = i;
1958                 }
1959             }
1960             if(j != imax){
1961                 for(int k=0;k<n;k++){
1962                     float dum = A[imax][k];
1963                     A[imax][k] = A[j][k];
1964                     A[j][k] = dum;
1965                 }
1966                 d = -d;
1967                 vv[imax] = vv[j];
1968             }
1969             index[j] = imax;
1970
1971             if(A[j][j] == 0.0){ A[j][j] = tiny; }
1972
1973             if(j != n-1){
1974                 float dum = 1.0/A[j][j];
1975                 for(int i=j+1;i<n;i++){ A[i][j] *= dum; }
1976             }
1977         }
1978     }
1979         catch(exception& e) {
1980                 m->errorOut(e, "LinearAlgebra", "ludcmp");
1981                 exit(1);
1982         }
1983 }
1984
1985 /*********************************************************************************************************************************/
1986
1987 void LinearAlgebra::lubksb(vector<vector<float> >& A, vector<int>& index, vector<float>& b){
1988     try {
1989         float total;
1990         int n = (int)A.size();
1991         int ii = 0;
1992
1993         for(int i=0;i<n;i++){
1994             if (m->control_pressed) { break; }
1995             int ip = index[i];
1996             total = b[ip];
1997             b[ip] = b[i];
1998
1999             if (ii != 0) {
2000                 for(int j=ii-1;j<i;j++){
2001                     total -= A[i][j] * b[j];
2002                 }
2003             }
2004             else if(total != 0){  ii = i+1;   }
2005             b[i] = total;
2006         }
2007         for(int i=n-1;i>=0;i--){
2008             total = b[i];
2009             for(int j=i+1;j<n;j++){ total -= A[i][j] * b[j];  }
2010             b[i] = total / A[i][i];
2011         }
2012     }
2013         catch(exception& e) {
2014                 m->errorOut(e, "LinearAlgebra", "lubksb");
2015                 exit(1);
2016         }
2017 }
2018
2019 /*********************************************************************************************************************************/
2020
2021 vector<vector<double> > LinearAlgebra::getInverse(vector<vector<double> > matrix){
2022     try {
2023         int n = (int)matrix.size();
2024
2025         vector<vector<double> > inverse(n);
2026         for(int i=0;i<n;i++){   inverse[i].assign(n, 0.0000);   }
2027
2028         vector<double> column(n, 0.0000);
2029         vector<int> index(n, 0);
2030         double dummy;
2031
2032         ludcmp(matrix, index, dummy);
2033
2034         for(int j=0;j<n;j++){
2035             if (m->control_pressed) { break; }
2036
2037             column.assign(n, 0);
2038
2039             column[j] = 1.0000;
2040
2041             lubksb(matrix, index, column);
2042
2043             for(int i=0;i<n;i++){   inverse[i][j] = column[i];  }
2044         }
2045
2046         return inverse;
2047     }
2048         catch(exception& e) {
2049                 m->errorOut(e, "LinearAlgebra", "getInverse");
2050                 exit(1);
2051         }
2052 }
2053 /*********************************************************************************************************************************/
2054 //modelled R lda function - MASS:::lda.default
2055 vector< vector<double> > LinearAlgebra::lda(vector< vector<double> >& a, vector<string> groups, vector< vector<double> >& means, bool& ignore) {
2056     try {
2057
2058         set<string> uniqueGroups;
2059         for (int i = 0; i < groups.size(); i++) { uniqueGroups.insert(groups[i]); }
2060         int numGroups = uniqueGroups.size();
2061
2062         map<string, int> quickIndex; //className to index. hoping to save counts, proportions and means in vectors to save time. This map will allow us to know index 0 in counts refers to group1.
2063         int count = 0;
2064         for (set<string>::iterator it = uniqueGroups.begin(); it != uniqueGroups.end(); it++) { quickIndex[*it] = count; count++; }
2065
2066         int numSampled = groups.size(); //number of sampled groups
2067         int numOtus = a.size(); //number of flagged bins
2068
2069         //counts <- as.vector(table(g)) //number of samples from each class in random sampling
2070         vector<int> counts; counts.resize(numGroups, 0);
2071         for (int i = 0; i < groups.size(); i++) {
2072             counts[quickIndex[groups[i]]]++;
2073         }
2074
2075         vector<double> proportions; proportions.resize(numGroups, 0.0);
2076         for (int i = 0; i < numGroups; i++) {  proportions[i] = counts[i] / (double) numSampled; }
2077
2078         means.clear(); //means[0] -> means[0][0] average for [group0][OTU0].
2079         means.resize(numGroups); for (int i = 0; i < means.size(); i++) { means[i].resize(numOtus, 0.0); }
2080         for (int j = 0; j < numSampled; j++) { //total for each class for each OTU
2081             for (int i = 0; i < numOtus; i++) { means[quickIndex[groups[j]]][i] += a[i][j]; }
2082         }
2083         //average for each class for each OTU
2084         for (int j = 0; j < numGroups; j++) { for (int i = 0; i < numOtus; i++) { means[j][i] /= counts[j]; }  }
2085
2086         //randCov <- x - group.means[g, ]
2087         vector< vector<double> > randCov; //randCov[0][0] -> (random sample value0 for OTU0 - average for samples group in OTU0). example OTU0, random sample 0.01 from class early. average of class early for OTU0 is 0.005. randCov[0][0] = (0.01-0.005)
2088         for (int i = 0; i < numOtus; i++) { //for each flagged OTU
2089             vector<double> tempRand;
2090             for (int j = 0; j < numSampled; j++) { tempRand.push_back(a[i][j] - means[quickIndex[groups[j]]][i]);  }
2091             randCov.push_back(tempRand);
2092         }
2093
2094         //find variance and std for each OTU
2095         //f1 <- sqrt(diag(var(x - group.means[g, ])))
2096         vector<double> stdF1;
2097         vector<double> ave;
2098         for (int i = 0; i < numOtus; i++) {
2099             stdF1.push_back(0.0);
2100             ave.push_back(m->getAverage(randCov[i]));
2101         }
2102
2103         for (int i = 0; i < numOtus; i++) {
2104             for (int j = 0; j < numSampled; j++) { stdF1[i] += ((randCov[i][j] - ave[i]) * (randCov[i][j] - ave[i]));  }
2105         }
2106
2107         //fac <- 1/(n - ng)
2108         double fac = 1 / (double) (numSampled-numGroups);
2109
2110         for (int i = 0; i < stdF1.size(); i++) {
2111             stdF1[i] /= (double) (numSampled-1);
2112             stdF1[i] = sqrt(stdF1[i]);
2113         }
2114
2115         vector< vector<double> > scaling; //[numOTUS][numOTUS]
2116         for (int i = 0; i < numOtus; i++) {
2117             vector<double> temp;
2118             for (int j = 0; j < numOtus; j++) {
2119                 if (i == j) { temp.push_back(1.0/stdF1[i]); }
2120                 else { temp.push_back(0.0); }
2121
2122             }
2123             scaling.push_back(temp);
2124         }
2125         /*
2126          cout << "scaling = " << endl;
2127          for (int i = 0; i < scaling.size(); i++) {
2128          for (int j = 0; j < scaling[i].size(); j++) { cout << scaling[i][j] << '\t'; }
2129          cout << endl;
2130          }*/
2131
2132         //X <- sqrt(fac) * ((x - group.means[g, ]) %*% scaling)
2133         vector< vector<double> > X = randCov; //[numOTUS][numSampled]
2134         //((x - group.means[g, ]) %*% scaling)
2135         //matrix multiplication of randCov and scaling
2136         LinearAlgebra linear;
2137         X = linear.matrix_mult(scaling, randCov); //[numOTUS][numOTUS] * [numOTUS][numSampled] = [numOTUS][numSampled]
2138         fac = sqrt(fac);
2139
2140         for (int i = 0; i < X.size(); i++) {
2141             for (int j = 0; j < X[i].size(); j++) { X[i][j] *= fac;  }
2142         }
2143
2144         vector<double> d;
2145         vector< vector<double> > v;
2146         vector< vector<double> > Xcopy; //X = [numOTUS][numSampled]
2147         bool transpose = false; //svd requires rows < columns, so if they are not then I need to transpose and look for the results in v.
2148         if (X.size() < X[0].size()) { Xcopy = linear.transpose(X); transpose=true; }
2149         else                        { Xcopy = X;                    }
2150         linear.svd(Xcopy, d, v); //Xcopy gets the results we want for v below, because R's version is [numSampled][numOTUS]
2151
2152         /*cout << "Xcopy = " << endl;
2153         for (int i = 0; i < Xcopy.size(); i++) {
2154             for (int j = 0; j < Xcopy[i].size(); j++) { cout << Xcopy[i][j] << '\t'; }
2155             cout << endl;
2156         }
2157         cout << "v = " << endl;
2158         for (int i = 0; i < v.size(); i++) {
2159             for (int j = 0; j < v[i].size(); j++) { cout << v[i][j] << '\t'; }
2160             cout << endl;
2161         }
2162          */
2163
2164         int rank = 0;
2165         set<int> goodColumns;
2166         //cout << "d = " << endl;
2167         for (int i = 0; i < d.size(); i++) {  if (d[i] > 0.0000000001) { rank++; goodColumns.insert(i); } } //cout << d[i] << endl;
2168
2169         if (rank == 0) {
2170             ignore=true; //m->mothurOut("[ERROR]: rank = 0: variables are numerically const\n"); m->control_pressed = true;
2171             return scaling; }
2172
2173         //scaling <- scaling %*% X.s$v[, 1L:rank] %*% diag(1/X.s$d[1L:rank], , rank)
2174         //X.s$v[, 1L:rank] = columns in Xcopy that correspond to "good" d values
2175         //diag(1/X.s$d[1L:rank], , rank) = matrix size rank * rank where the diagonal is 1/"good" dvalues
2176         /*example:
2177          d
2178          [1] 3.721545e+00 3.034607e+00 2.296649e+00 7.986927e-16 6.922408e-16
2179          [6] 5.471102e-16
2180
2181          $v
2182          [,1]        [,2]        [,3]        [,4]        [,5]        [,6]
2183          [1,]  0.31122175  0.10944725  0.20183340 -0.30136820  0.60786235 -0.13537095
2184          [2,] -0.29563726 -0.20568893  0.11233366 -0.05073289  0.48234270  0.21965978
2185          ...
2186
2187          [1] "X.s$v[, 1L:rank]"
2188          [,1]        [,2]        [,3]
2189          [1,]  0.31122175  0.10944725  0.20183340
2190          [2,] -0.29563726 -0.20568893  0.11233366
2191          ...
2192          [1] "1/X.s$d[1L:rank]"
2193          [1] 0.2687056 0.3295320 0.4354170
2194
2195          [1] "diag(1/X.s$d[1L:rank], , rank)"
2196          [,1]     [,2]     [,3]
2197          [1,] 0.2687056 0.000000 0.000000
2198          [2,] 0.0000000 0.329532 0.000000
2199          [3,] 0.0000000 0.000000 0.435417
2200          */
2201         if (transpose) {
2202             Xcopy = linear.transpose(v);
2203             /*
2204             cout << "Xcopy = " << endl;
2205             for (int i = 0; i < Xcopy.size(); i++) {
2206                 for (int j = 0; j < Xcopy[i].size(); j++) { cout << Xcopy[i][j] << '\t'; }
2207                 cout << endl;
2208             }*/
2209         }
2210         v.clear(); //store "good" columns - X.s$v[, 1L:rank]
2211         v.resize(Xcopy.size()); //[numOTUS]["good" columns]
2212         for (set<int>::iterator it = goodColumns.begin(); it != goodColumns.end(); it++) {
2213             for (int i = 0; i < Xcopy.size(); i++) {
2214                 v[i].push_back(Xcopy[i][*it]);
2215             }
2216         }
2217
2218         vector< vector<double> > diagRanks; diagRanks.resize(rank);
2219         for (int i = 0; i < rank; i++) { diagRanks[i].resize(rank, 0.0); }
2220         count = 0;
2221         for (set<int>::iterator it = goodColumns.begin(); it != goodColumns.end(); it++) {  diagRanks[count][count] = 1.0 / d[*it]; count++; }
2222
2223         scaling = linear.matrix_mult(linear.matrix_mult(scaling, v), diagRanks); //([numOTUS][numOTUS]*[numOTUS]["good" columns]) = [numOTUS]["good" columns] then ([numOTUS]["good" columns] * ["good" columns]["good" columns] = scaling = [numOTUS]["good" columns]
2224
2225         /*cout << "scaling = " << endl;
2226         for (int i = 0; i < scaling.size(); i++) {
2227             for (int j = 0; j < scaling[i].size(); j++) { cout << scaling[i][j] << '\t'; }
2228             cout << endl;
2229         }*/
2230
2231         //Note: linear.matrix_mult [1][numGroups] * [numGroups][numOTUs] - columns in first must match rows in second, returns matrix[1][numOTUs]
2232         vector< vector<double> > prior; prior.push_back(proportions);
2233         vector< vector<double> >  xbar = linear.matrix_mult(prior, means);
2234         vector<double> xBar = xbar[0]; //length numOTUs
2235
2236         /*cout << "xbar" << endl;
2237         for (int j = 0; j < numOtus; j++) {  cout << xBar[j] <<'\t'; }cout <<  endl;*/
2238         //fac <- 1/(ng - 1)
2239         fac = 1 / (double) (numGroups-1);
2240         //scale(group.means, center = xbar, scale = FALSE) %*% scaling
2241         vector< vector<double> > scaledMeans = means; //[numGroups][numOTUs]
2242         for (int i = 0; i < numGroups; i++) {
2243             for (int j = 0; j < numOtus; j++) {  scaledMeans[i][j] -= xBar[j]; }
2244         }
2245         scaledMeans = linear.matrix_mult(scaledMeans, scaling); //[numGroups][numOTUS]*[numOTUS]["good"columns] = [numGroups]["good"columns]
2246
2247
2248         //sqrt((n * prior) * fac)
2249         vector<double> temp = proportions; //[numGroups]
2250         for (int i = 0; i < temp.size(); i++) { temp[i] *= numSampled * fac; temp[i] = sqrt(temp[i]);  }
2251
2252         //X <- sqrt((n * prior) * fac) * (scale(group.means, center = xbar, scale = FALSE) %*% scaling)
2253         //X <- temp * scaledMeans
2254         X.clear(); X = scaledMeans; //[numGroups]["good"columns]
2255         for (int i = 0; i < X.size(); i++) {
2256             for (int j = 0; j < X[i].size(); j++) {  X[i][j] *= temp[j];  }
2257         }
2258         /*
2259         cout << "X = " << endl;
2260         for (int i = 0; i < X.size(); i++) {
2261             for (int j = 0; j < X[i].size(); j++) { cout << X[i][j] << '\t'; }
2262             cout << endl;
2263         }
2264         */
2265
2266         d.clear(); v.clear();
2267         //we want to transpose so results are in Xcopy, but if that makes rows > columns then we don't since svd requires rows < cols.
2268         transpose=false;
2269         if (X.size() > X[0].size()) {   Xcopy = X;  transpose=true;     }
2270         else                        {   Xcopy = linear.transpose(X);    }
2271         linear.svd(Xcopy, d, v); //Xcopy gets the results we want for v below
2272         /*cout << "Xcopy = " << endl;
2273         for (int i = 0; i < Xcopy.size(); i++) {
2274             for (int j = 0; j < Xcopy[i].size(); j++) { cout << Xcopy[i][j] << '\t'; }
2275             cout << endl;
2276         }
2277
2278         cout << "v = " << endl;
2279         for (int i = 0; i < v.size(); i++) {
2280             for (int j = 0; j < v[i].size(); j++) { cout << v[i][j] << '\t'; }
2281             cout << endl;
2282         }
2283
2284         cout << "d = " << endl;
2285         for (int i = 0; i < d.size(); i++) { cout << d[i] << endl; }*/
2286
2287         //rank <- sum(X.s$d > tol * X.s$d[1L])
2288         //X.s$d[1L] = larger value in d vector
2289         double largeD = m->max(d);
2290         rank = 0; goodColumns.clear();
2291         for (int i = 0; i < d.size(); i++) { if (d[i] > (0.0000000001*largeD)) { rank++; goodColumns.insert(i); } }
2292
2293         if (rank == 0) {
2294             ignore=true;//m->mothurOut("[ERROR]: rank = 0: class means are numerically identical.\n"); m->control_pressed = true;
2295             return scaling; }
2296
2297         if (transpose) { Xcopy = linear.transpose(v);  }
2298         //scaling <- scaling %*% X.s$v[, 1L:rank] - scaling * "good" columns
2299         v.clear(); //store "good" columns - X.s$v[, 1L:rank]
2300         v.resize(Xcopy.size()); //Xcopy = ["good"columns][numGroups]
2301         for (set<int>::iterator it = goodColumns.begin(); it != goodColumns.end(); it++) {
2302             for (int i = 0; i < Xcopy.size(); i++) {
2303                 v[i].push_back(Xcopy[i][*it]);
2304             }
2305         }
2306
2307         scaling = linear.matrix_mult(scaling, v); //[numOTUS]["good" columns] * ["good"columns][new "good" columns]
2308
2309         /*cout << "scaling = " << endl;
2310         for (int i = 0; i < scaling.size(); i++) {
2311             for (int j = 0; j < scaling[i].size(); j++) { cout << scaling[i][j] << '\t'; }
2312             cout << endl;
2313         }*/
2314         ignore=false;
2315         return scaling;
2316     }
2317         catch(exception& e) {
2318                 m->errorOut(e, "LinearAlgebra", "lda");
2319                 exit(1);
2320         }
2321 }
2322 /*********************************************************************************************************************************/
2323 //Singular value decomposition (SVD) - adapted from http://svn.lirec.eu/libs/magicsquares/src/SVD.cpp
2324 /*
2325  * svdcomp - SVD decomposition routine.
2326  * Takes an mxn matrix a and decomposes it into udv, where u,v are
2327  * left and right orthogonal transformation matrices, and d is a
2328  * diagonal matrix of singular values.
2329  *
2330  * This routine is adapted from svdecomp.c in XLISP-STAT 2.1 which is
2331  * code from Numerical Recipes adapted by Luke Tierney and David Betz.
2332  *
2333  * Input to dsvd is as follows:
2334  *   a = mxn matrix to be decomposed, gets overwritten with u
2335  *   m = row dimension of a
2336  *   n = column dimension of a
2337  *   w = returns the vector of singular values of a
2338  *   v = returns the right orthogonal transformation matrix
2339  */
2340
2341 int LinearAlgebra::svd(vector< vector<double> >& a, vector<double>& w, vector< vector<double> >& v) {
2342     try {
2343         int flag, i, its, j, jj, k, l, nm;
2344         double c, f, h, s, x, y, z;
2345         double anorm = 0.0, g = 0.0, scale = 0.0;
2346
2347         int numRows = a.size(); if (numRows == 0) { return 0; }
2348         int numCols = a[0].size();
2349         w.resize(numCols, 0.0);
2350         v.resize(numCols); for (int i = 0; i < numCols; i++) { v[i].resize(numRows, 0.0); }
2351
2352         vector<double> rv1; rv1.resize(numCols, 0.0);
2353         if (numRows < numCols){  m->mothurOut("[ERROR]: numRows < numCols\n"); m->control_pressed = true; return 0; }
2354
2355         /* Householder reduction to bidiagonal form */
2356         for (i = 0; i < numCols; i++)
2357         {
2358             /* left-hand reduction */
2359             l = i + 1;
2360             rv1[i] = scale * g;
2361             g = s = scale = 0.0;
2362             if (i < numRows)
2363             {
2364                 for (k = i; k < numRows; k++)
2365                     scale += fabs((double)a[k][i]);
2366                 if (scale)
2367                 {
2368                     for (k = i; k < numRows; k++)
2369                     {
2370                         a[k][i] = (double)((double)a[k][i]/scale);
2371                         s += ((double)a[k][i] * (double)a[k][i]);
2372                     }
2373                     f = (double)a[i][i];
2374                     g = -SIGN(sqrt(s), f);
2375                     h = f * g - s;
2376                     a[i][i] = (double)(f - g);
2377                     if (i != numCols - 1)
2378                     {
2379                         for (j = l; j < numCols; j++)
2380                         {
2381                             for (s = 0.0, k = i; k < numRows; k++)
2382                                 s += ((double)a[k][i] * (double)a[k][j]);
2383                             f = s / h;
2384                             for (k = i; k < numRows; k++)
2385                                 a[k][j] += (double)(f * (double)a[k][i]);
2386                         }
2387                     }
2388                     for (k = i; k < numRows; k++)
2389                         a[k][i] = (double)((double)a[k][i]*scale);
2390                 }
2391             }
2392             w[i] = (double)(scale * g);
2393
2394             /* right-hand reduction */
2395             g = s = scale = 0.0;
2396             if (i < numRows && i != numCols - 1)
2397             {
2398                 for (k = l; k < numCols; k++)
2399                     scale += fabs((double)a[i][k]);
2400                 if (scale)
2401                 {
2402                     for (k = l; k < numCols; k++)
2403                     {
2404                         a[i][k] = (double)((double)a[i][k]/scale);
2405                         s += ((double)a[i][k] * (double)a[i][k]);
2406                     }
2407                     f = (double)a[i][l];
2408                     g = -SIGN(sqrt(s), f);
2409                     h = f * g - s;
2410                     a[i][l] = (double)(f - g);
2411                     for (k = l; k < numCols; k++)
2412                         rv1[k] = (double)a[i][k] / h;
2413                     if (i != numRows - 1)
2414                     {
2415                         for (j = l; j < numRows; j++)
2416                         {
2417                             for (s = 0.0, k = l; k < numCols; k++)
2418                                 s += ((double)a[j][k] * (double)a[i][k]);
2419                             for (k = l; k < numCols; k++)
2420                                 a[j][k] += (double)(s * rv1[k]);
2421                         }
2422                     }
2423                     for (k = l; k < numCols; k++)
2424                         a[i][k] = (double)((double)a[i][k]*scale);
2425                 }
2426             }
2427             anorm = max(anorm, (fabs((double)w[i]) + fabs(rv1[i])));
2428         }
2429
2430         /* accumulate the right-hand transformation */
2431         for (i = numCols - 1; i >= 0; i--)
2432         {
2433             if (i < numCols - 1)
2434             {
2435                 if (g)
2436                 {
2437                     for (j = l; j < numCols; j++)
2438                         v[j][i] = (double)(((double)a[i][j] / (double)a[i][l]) / g);
2439                     /* double division to avoid underflow */
2440                     for (j = l; j < numCols; j++)
2441                     {
2442                         for (s = 0.0, k = l; k < numCols; k++)
2443                             s += ((double)a[i][k] * (double)v[k][j]);
2444                         for (k = l; k < numCols; k++)
2445                             v[k][j] += (double)(s * (double)v[k][i]);
2446                     }
2447                 }
2448                 for (j = l; j < numCols; j++)
2449                     v[i][j] = v[j][i] = 0.0;
2450             }
2451             v[i][i] = 1.0;
2452             g = rv1[i];
2453             l = i;
2454         }
2455
2456         /* accumulate the left-hand transformation */
2457         for (i = numCols - 1; i >= 0; i--)
2458         {
2459             l = i + 1;
2460             g = (double)w[i];
2461             if (i < numCols - 1)
2462                 for (j = l; j < numCols; j++)
2463                     a[i][j] = 0.0;
2464             if (g)
2465             {
2466                 g = 1.0 / g;
2467                 if (i != numCols - 1)
2468                 {
2469                     for (j = l; j < numCols; j++)
2470                     {
2471                         for (s = 0.0, k = l; k < numRows; k++)
2472                             s += ((double)a[k][i] * (double)a[k][j]);
2473                         f = (s / (double)a[i][i]) * g;
2474                         for (k = i; k < numRows; k++)
2475                             a[k][j] += (double)(f * (double)a[k][i]);
2476                     }
2477                 }
2478                 for (j = i; j < numRows; j++)
2479                     a[j][i] = (double)((double)a[j][i]*g);
2480             }
2481             else
2482             {
2483                 for (j = i; j < numRows; j++)
2484                     a[j][i] = 0.0;
2485             }
2486             ++a[i][i];
2487         }
2488
2489         /* diagonalize the bidiagonal form */
2490         for (k = numCols - 1; k >= 0; k--)
2491         {                             /* loop over singular values */
2492             for (its = 0; its < 30; its++)
2493             {                         /* loop over allowed iterations */
2494                 flag = 1;
2495                 for (l = k; l >= 0; l--)
2496                 {                     /* test for splitting */
2497                     nm = l - 1;
2498                     if (fabs(rv1[l]) + anorm == anorm)
2499                     {
2500                         flag = 0;
2501                         break;
2502                     }
2503                     if (fabs((double)w[nm]) + anorm == anorm)
2504                         break;
2505                 }
2506                 if (flag)
2507                 {
2508                     c = 0.0;
2509                     s = 1.0;
2510                     for (i = l; i <= k; i++)
2511                     {
2512                         f = s * rv1[i];
2513                         if (fabs(f) + anorm != anorm)
2514                         {
2515                             g = (double)w[i];
2516                             h = pythag(f, g);
2517                             w[i] = (double)h;
2518                             h = 1.0 / h;
2519                             c = g * h;
2520                             s = (- f * h);
2521                             for (j = 0; j < numRows; j++)
2522                             {
2523                                 y = (double)a[j][nm];
2524                                 z = (double)a[j][i];
2525                                 a[j][nm] = (double)(y * c + z * s);
2526                                 a[j][i] = (double)(z * c - y * s);
2527                             }
2528                         }
2529                     }
2530                 }
2531                 z = (double)w[k];
2532                 if (l == k)
2533                 {                  /* convergence */
2534                     if (z < 0.0)
2535                     {              /* make singular value nonnegative */
2536                         w[k] = (double)(-z);
2537                         for (j = 0; j < numCols; j++)
2538                             v[j][k] = (-v[j][k]);
2539                     }
2540                     break;
2541                 }
2542                 if (its >= 30) {
2543                     m->mothurOut("No convergence after 30,000! iterations \n"); m->control_pressed = true;
2544                     return(0);
2545                 }
2546
2547                 /* shift from bottom 2 x 2 minor */
2548                 x = (double)w[l];
2549                 nm = k - 1;
2550                 y = (double)w[nm];
2551                 g = rv1[nm];
2552                 h = rv1[k];
2553                 f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
2554                 g = pythag(f, 1.0);
2555                 f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x;
2556
2557                 /* next QR transformation */
2558                 c = s = 1.0;
2559                 for (j = l; j <= nm; j++)
2560                 {
2561                     i = j + 1;
2562                     g = rv1[i];
2563                     y = (double)w[i];
2564                     h = s * g;
2565                     g = c * g;
2566                     z = pythag(f, h);
2567                     rv1[j] = z;
2568                     c = f / z;
2569                     s = h / z;
2570                     f = x * c + g * s;
2571                     g = g * c - x * s;
2572                     h = y * s;
2573                     y = y * c;
2574                     for (jj = 0; jj < numCols; jj++)
2575                     {
2576                         x = (double)v[jj][j];
2577                         z = (double)v[jj][i];
2578                         v[jj][j] = (float)(x * c + z * s);
2579                         v[jj][i] = (float)(z * c - x * s);
2580                     }
2581                     z = pythag(f, h);
2582                     w[j] = (float)z;
2583                     if (z)
2584                     {
2585                         z = 1.0 / z;
2586                         c = f * z;
2587                         s = h * z;
2588                     }
2589                     f = (c * g) + (s * y);
2590                     x = (c * y) - (s * g);
2591                     for (jj = 0; jj < numRows; jj++)
2592                     {
2593                         y = (double)a[jj][j];
2594                         z = (double)a[jj][i];
2595                         a[jj][j] = (double)(y * c + z * s);
2596                         a[jj][i] = (double)(z * c - y * s);
2597                     }
2598                 }
2599                 rv1[l] = 0.0;
2600                 rv1[k] = f;
2601                 w[k] = (double)x;
2602             }
2603         }
2604
2605         return(0);
2606
2607     }
2608         catch(exception& e) {
2609                 m->errorOut(e, "LinearAlgebra", "svd");
2610                 exit(1);
2611         }
2612 }
2613
2614