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