mothurmetastats.cpp

   1 /*
   2  *  mothurmetastats.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 7/6/11.
   6  *  Copyright 2011 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "mothurmetastats.h"
  11 #include "mothurfisher.h"
  12
  13 /***********************************************************/
  14 MothurMetastats::MothurMetastats(double t, int n) {
  15         try {
  16                 m = MothurOut::getInstance();
  17                 threshold = t;
  18                 numPermutations = n;
  19
  20         }catch(exception& e) {
  21                 m->errorOut(e, "MothurMetastats", "MothurMetastats");
  22                 exit(1);
  23         }
  24 }
  25 /***********************************************************/
  26 MothurMetastats::~MothurMetastats() {}
  27 /***********************************************************/
  28 //main metastats function
  29 int MothurMetastats::runMetastats(string outputFileName, vector< vector<double> >& data, int secondGroupingStart) {
  30         try {
  31                 int bflag = 0;
  32                 row = data.size();               //numBins
  33                 column = data[0].size(); //numGroups in subset
  34                 int size = row*column;
  35
  36                 //consistent with original, but this should never be true
  37                 if ((secondGroupingStart >= column) || (secondGroupingStart <= 0)) { m->mothurOut("[ERROR]: Check your g value."); m->mothurOutEndLine(); return 0; }
  38
  39                 //Initialize the matrices
  40                 vector<double> pmatrix; pmatrix.resize(size, 0.0);
  41                 vector<double> permuted; permuted.resize(size, 0.0);
  42                 vector< vector<double> > storage; storage.resize(row);
  43                 for (int i = 0; i < storage.size(); i++) { storage[i].resize(9, 0.0); }
  44
  45                 //Produces the sum of each column
  46                 vector<double> total; total.resize(column, 0.0);
  47                 vector<double> ratio; ratio.resize(column, 0.0);
  48                 double total1 = 0.0; double total2 = 0.0;
  49
  50                 //total[i] = total abundance for group[i]
  51                 for (int i = 0; i < column; i++) {
  52                         for (int j = 0; j < row; j++) {
  53                                 total[i] += data[j][i];
  54                         }
  55                 }
  56
  57                 //total for first grouping
  58                 for (int i = 0; i < secondGroupingStart; i++) { total1 += total[i]; }
  59
  60                 //total for second grouping
  61                 for (int i = secondGroupingStart; i < column; i++) { total2 += total[i]; }
  62
  63                 //Creates the ratios by first finding the minimum of totals
  64                 double min = total[0];
  65                 for (int i = 0; i < total.size(); i++) {
  66                          if (total[i] < min) { min = total[i]; }
  67                 }
  68
  69                 //sanity check
  70                 if (min <= 0.0) { m->mothurOut("[ERROR]: the sum of one of the columns <= 0."); m->mothurOutEndLine(); return 0; }
  71
  72                 //Ratio time...
  73                 for(int i = 0; i < ratio.size(); i++){  ratio[i] = total[i] / min; }
  74
  75                 //Change matrix into an array as received by R for compatibility - kept to be consistent with original
  76                 int count = 0;
  77                 for(int i = 0; i < column; i++){
  78                         for(int j = 0; j < row; j++){
  79                                 pmatrix[count]=data[j][i];
  80                                 count++;
  81                         }
  82                 }
  83
  84                 if(row == 1){
  85                         for (int i =0; i < column; i++){ pmatrix[i] /= ratio[i]; }
  86                 }else {
  87                         count = 0; int j=-1;
  88
  89                         for (int i=0; i < size; i++) {
  90                                 if (count % row == 0) { j++; }
  91                                 pmatrix[i] /= ratio[j];
  92                                 count++;
  93                         }
  94                 }
  95
  96                 vector<double> permuted_ttests; permuted_ttests.resize(row, 0.0);
  97                 vector<double> pvalues;                 pvalues.resize(row, 0.0);
  98                 vector<double> tinitial;                tinitial.resize(row, 0.0);
  99
 100                 if (m->control_pressed) { return 1; }
 101
 102                 //Find the initial values for the matrix.
 103                 start(pmatrix, secondGroupingStart, tinitial, storage);
 104
 105                 if (m->control_pressed) { return 1; }
 106
 107                 // Start the calculations.
 108                 if ( (column == 2) || (secondGroupingStart < 8) || ((column-secondGroupingStart) < 8) ){
 109
 110                         vector<double> fish;    fish.resize(row, 0.0);
 111                         vector<double> fish2;   fish2.resize(row, 0.0);
 112
 113                         for(int i = 0; i < row; i++){
 114
 115                                 for(int j = 0; j < secondGroupingStart; j++)            { fish[i] += data[i][j];        }
 116                                 for(int j = secondGroupingStart; j < column; j++)       { fish2[i] += data[i][j];       }
 117
 118                                 //vector<double> tempData; tempData.resize(4, 0.0);
 119                                 double f11, f12, f21, f22;
 120                                 f11 = fish[i];
 121                                 f12 = fish2[i];
 122                                 f21 = total1 - fish[i];
 123                                 f22 = total2 - fish2[i];
 124
 125                                 double pre = 0.0;
 126
 127                                 MothurFisher fisher;
 128                                 pre = fisher.fexact(f11, f12, f21, f22);
 129
 130                                 if (m->control_pressed) { return 1; }
 131
 132                                 if (pre > 0.999999999)  { pre = 1.0; }
 133                                 storage[i][8] = pre;
 134                                 pvalues[i] = pre;
 135                         }
 136
 137                 }else {
 138
 139                         testp(permuted_ttests, permuted, pmatrix, secondGroupingStart, tinitial, pvalues);
 140
 141                         if (m->control_pressed) { return 1; }
 142
 143                         // Checks to make sure the matrix isn't sparse.
 144                         vector<double> sparse;          sparse.resize(row, 0.0);
 145                         vector<double> sparse2;         sparse2.resize(row, 0.0);
 146
 147                         int c = 0;
 148
 149                         for(int i = 0; i < row; i++){
 150
 151                                 for(int j = 0; j < secondGroupingStart; j++)    {       sparse[i] += data[i][j];        }
 152                                 if(sparse[i] < (double)secondGroupingStart)             {       c++;                                            }
 153
 154                                 // ?<= for col
 155                                 for(int j = secondGroupingStart; j < column; j++)               {  sparse2[i] += data[i][j]; }
 156                                 if( (sparse2[i] < (double)(column-secondGroupingStart)))        { c++;                                           }
 157
 158                                 if (c == 2) {
 159                                         c=0;
 160                                         double f11,f12,f21,f22;
 161
 162                                         f11=sparse[i];  sparse[i]=0;
 163                                         f12=sparse2[i];  sparse2[i]=0;
 164                                         f21 = total1 - f11;
 165                                         f22 = total2 - f12;
 166
 167                                         double pre = 0.0;
 168
 169                                         MothurFisher fisher;
 170                                         pre = fisher.fexact(f11, f12, f21, f22);
 171
 172                                         if (m->control_pressed) { return 1; }
 173
 174                                         if (pre > 0.999999999){
 175                                                 pre = 1.0;
 176                                         }
 177
 178                                         storage[i][8] = pre;
 179                                         pvalues[i] = pre;
 180                                 }
 181                         }
 182
 183                         bflag = 1;
 184                 }
 185
 186                 // Calculates the mean of counts (not normalized)
 187                 vector< vector<double> > temp; temp.resize(row);
 188                 for (int i = 0; i < temp.size(); i++) { temp[i].resize(2, 0.0); }
 189
 190                 for (int j = 0; j < row; j++){
 191                         if (m->control_pressed) { return 1; }
 192
 193                         for (int i = 0; i < secondGroupingStart; i++){ temp[j][0] += data[j][i]; }
 194                         temp[j][0] /= (double)secondGroupingStart;
 195
 196                         for(int i = secondGroupingStart; i < column; i++){ temp[j][1] += data[j][i]; }
 197                         temp[j][1] /= (double)(column-secondGroupingStart);
 198                 }
 199
 200                 for(int i = 0; i < row; i++){
 201                         if (m->control_pressed) { return 1; }
 202
 203                         storage[i][3]=temp[i][0];
 204                         storage[i][7]=temp[i][1];
 205                         storage[i][8]=pvalues[i];
 206                 }
 207                 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
 208                 cout << "pvalues" << endl;
 209                 for (int i = 0; i < row; i++){ if (pvalues[i] < 0.0000001) {
 210                         pvalues[i] = 0.0;
 211                 }cout << pvalues[i] << '\t'; }
 212                 cout << endl;
 213                 calc_qvalues(pvalues);
 214
 215                 // BACKUP checks
 216                 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
 217                 for (int i = 0; i < row; i++){
 218
 219                         if (m->control_pressed) { return 1; }
 220
 221                         if(pvalues[i] < threshold){
 222                                 m->mothurOut("Feature " + toString((i+1)) + " is significant, p = ");
 223                                 cout << pvalues[i];
 224                                 m->mothurOutJustToLog(toString(pvalues[i])); m->mothurOutEndLine();
 225                         }
 226                 }
 227
 228                 // And now we write the files to a text file.
 229                 struct tm *local;
 230                 time_t t; t = time(NULL);
 231                 local = localtime(&t);
 232
 233                 ofstream out;
 234                 m->openOutputFile(outputFileName, out);
 235                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 236
 237                 out << "Local time and date of test: " << asctime(local) << endl;
 238                 out << "# rows = " << row << ", # col = " << column << ", g = " << secondGroupingStart << endl << endl;
 239                 if (bflag == 1){ out << numPermutations << " permutations" << endl << endl;     }
 240
 241                 //output column headings - not really sure... documentation labels 9 columns, there are 10 in the output file
 242                 //storage 0 = meanGroup1 - line 529, 1 = varGroup1 - line 532, 2 = err rate1 - line 534, 3 = mean of counts group1?? - line 291, 4 = meanGroup2 - line 536, 5 = varGroup2 - line 539, 6 = err rate2 - line 541, 7 = mean of counts group2?? - line 292, 8 = pvalues - line 293
 243                 out << "OTU\tmean(group1)\tvariance(group1)\tstderr(group1)\tmean_of_counts(group1)\tmean(group2)\tvariance(group2)\tstderr(group2)\tmean_of_counts(group1)\tp-value\n";
 244
 245                 for(int i = 0; i < row; i++){
 246                         if (m->control_pressed) { out.close(); return 0; }
 247
 248                         out << (i+1);
 249                         for(int j = 0; j < 9; j++){ out << '\t' << storage[i][j]; }
 250                         out << endl;
 251                 }
 252
 253                 out << endl << endl;
 254                 out.close();
 255
 256
 257                 return 0;
 258
 259         }catch(exception& e) {
 260                 m->errorOut(e, "MothurMetastats", "runMetastats");
 261                 exit(1);
 262         }
 263 }
 264 /***********************************************************/
 265 //Find the initial values for the matrix
 266 int MothurMetastats::start(vector<double>& Imatrix, int secondGroupingStart, vector<double>& initial, vector< vector<double> >& storage) {
 267         try {
 268
 269                 int a = row; a*=4;
 270
 271                 double xbardiff = 0.0; double denom = 0.0;
 272                 vector<double> store;   store.resize(a, 0.0);
 273                 vector<double> tool;    tool.resize(a, 0.0);
 274                 vector< vector<double> > C1; C1.resize(row);
 275                 for (int i = 0; i < C1.size(); i++) { C1[i].resize(3, 0.0); }
 276                 vector< vector<double> > C2; C2.resize(row);
 277                 for (int i = 0; i < C2.size(); i++) { C2[i].resize(3, 0.0); }
 278
 279                 meanvar(Imatrix, secondGroupingStart, store);
 280
 281                 if (m->control_pressed) { return 0; }
 282
 283                 //copy store into tool
 284                 tool = store;
 285
 286                 for (int i = 0; i < row; i++){
 287                         C1[i][0]=tool[i]; //mean group 1
 288                         storage[i][0]=C1[i][0];
 289                         C1[i][1]=tool[i+row+row]; // var group 1
 290                         storage[i][1]=C1[i][1];
 291                         C1[i][2]=C1[i][1]/(secondGroupingStart);
 292                         storage[i][2]=sqrt(C1[i][2]);
 293
 294                         C2[i][0]=tool[i+row]; // mean group 2
 295                         storage[i][4]=C2[i][0];
 296                         C2[i][1]=tool[i+row+row+row]; // var group 2
 297                         storage[i][5]=C2[i][1];
 298                         C2[i][2]=C2[i][1]/(column-secondGroupingStart);
 299                         storage[i][6]=sqrt(C2[i][2]);
 300                 }
 301
 302                 if (m->control_pressed) { return 0; }
 303
 304                 for (int i = 0; i < row; i++){
 305                         xbardiff = C1[i][0]-C2[i][0];
 306                         denom = sqrt(C1[i][2]+C2[i][2]);
 307                         initial[i]=fabs(xbardiff/denom);
 308                 }
 309
 310                 return 0;
 311
 312         }catch(exception& e) {
 313                 m->errorOut(e, "MothurMetastats", "start");
 314                 exit(1);
 315         }
 316 }
 317 /***********************************************************/
 318 int MothurMetastats::meanvar(vector<double>& pmatrix, int secondGroupingStart, vector<double>& store) {
 319         try {
 320                 vector<double> temp;    temp.resize(row, 0.0);
 321                 vector<double> temp2;   temp2.resize(row, 0.0);
 322                 vector<double> var;             var.resize(row, 0.0);
 323                 vector<double> var2;    var2.resize(row, 0.0);
 324
 325                 double a = secondGroupingStart;
 326                 double b = column - a;
 327                 int m = a * row;
 328                 int n = row * column;
 329
 330                 for (int i = 0; i < m; i++)             { temp[i%row] += pmatrix[i];    }
 331                 for (int i = 0; i < n; i++)             { temp2[i%row]+= pmatrix[i];    }
 332                 for (int i = 0; i < row; i++)   { temp2[i] -= temp[i];          }
 333                 for (int i = 0; i <= row-1;i++) {
 334                         store[i] = temp[i]/a;
 335                         store[i+row]=temp2[i]/b;
 336                 }
 337
 338                 //That completes the mean calculations.
 339
 340                 for (int i = 0; i < m; i++)             { var[i%row] += pow((pmatrix[i]-store[i%row]),2);               }
 341                 for (int i = m; i < n; i++)             { var2[i%row]+= pow((pmatrix[i]-store[(i%row)+row]),2); }
 342                 for (int i = 0; i <= row-1; i++){
 343                         store[i+2*row]=var[i]/(a-1);
 344                         store[i+3*row]=var2[i]/(b-1);
 345                 }
 346
 347                 // That completes var calculations.
 348
 349                 return 0;
 350
 351         }catch(exception& e) {
 352                 m->errorOut(e, "MothurMetastats", "meanvar");
 353                 exit(1);
 354         }
 355 }
 356 /***********************************************************/
 357 int MothurMetastats::testp(vector<double>& permuted_ttests, vector<double>& permuted, vector<double>& Imatrix, int secondGroupingStart, vector<double>& Tinitial, vector<double>& ps) {
 358         try {
 359
 360                 vector<double> Tvalues;         Tvalues.resize(row, 0.0);
 361                 vector<double> counter;         counter.resize(row, 0.0);
 362                 int a, b, n;
 363
 364                 a = numPermutations;
 365                 b = row;
 366                 n = a*b;
 367
 368                 for (int j = 1; j <= row; j++)  {
 369                         if (m->control_pressed) { return 0; }
 370                         permute_matrix(Imatrix, permuted, secondGroupingStart, Tvalues, Tinitial, counter);
 371                 }
 372
 373                 for(int j = 0; j < row; j++)    {
 374                         if (m->control_pressed) { return 0; }
 375                         ps[j] = ((counter[j]+1)/(double)(a+1));
 376                 }
 377
 378                 return 0;
 379
 380         }catch(exception& e) {
 381                 m->errorOut(e, "MothurMetastats", "testp");
 382                 exit(1);
 383         }
 384 }
 385 /***********************************************************/
 386 int MothurMetastats::permute_matrix(vector<double>& Imatrix, vector<double>& permuted, int secondGroupingStart, vector<double>& trial_ts, vector<double>& Tinitial, vector<double>& counter1){
 387         try {
 388
 389                 vector<int> y; y.resize(column, 0);
 390                 for (int i = 1; i <= column; i++){ y[i-1] = i; }
 391
 392                 permute_array(y);
 393
 394                 int f = 0; int c = 0; int k = 0;
 395                 for (int i = 0; i < column; i++){
 396
 397                         if (m->control_pressed) { return 0; }
 398
 399                         f = y[i]; //column number
 400                         c = 1;
 401                         c *= (f-1);
 402                         c *= row;
 403                         if (f == 1){ c = 0; } // starting value position in the Imatrix
 404
 405                         for(int j = 1; j <= row; j++){
 406                                 permuted[k] = Imatrix[c];
 407                                 c++; k++;
 408                         }
 409                 }
 410
 411                 calc_twosample_ts(permuted, secondGroupingStart, trial_ts, Tinitial, counter1);
 412
 413                 return 0;
 414
 415         }catch(exception& e) {
 416                 m->errorOut(e, "MothurMetastats", "permute_matrix");
 417                 exit(1);
 418         }
 419 }
 420 /***********************************************************/
 421 int MothurMetastats::permute_array(vector<int>& array) {
 422         try {
 423                 static int seeded = 0;
 424
 425                 if (! seeded) {
 426                         seeded = 1;
 427                         srand(time(NULL));
 428                 }
 429
 430                 for (int i = 0; i < array.size(); i++) {
 431                         if (m->control_pressed) { return 0; }
 432
 433                         int selection = rand() % (array.size() - i);
 434                         int tmp = array[i + selection];
 435                         array[i + selection] = array[i];
 436                         array[i] = tmp;
 437                 }
 438
 439                 return 0;
 440
 441         }catch(exception& e) {
 442                 m->errorOut(e, "MothurMetastats", "permute_array");
 443                 exit(1);
 444         }
 445 }
 446 /***********************************************************/
 447 int MothurMetastats::calc_twosample_ts(vector<double>& Pmatrix, int secondGroupingStart, vector<double>& Ts, vector<double>& Tinitial, vector<double>& counter) {
 448         try {
 449                 int a = row * 4;
 450
 451                 vector< vector<double> > C1; C1.resize(row);
 452                 for (int i = 0; i < C1.size(); i++) { C1[i].resize(3, 0.0); }
 453                 vector< vector<double> > C2; C2.resize(row);
 454                 for (int i = 0; i < C2.size(); i++) { C2[i].resize(3, 0.0); }
 455                 vector<double> storage; storage.resize(a, 0.0);
 456                 vector<double> tool;    tool.resize(a, 0.0);
 457                 double xbardiff = 0.0; double denom = 0.0;
 458
 459                 meanvar(Pmatrix, secondGroupingStart, storage);
 460
 461                 for(int i = 0;i <= (a-1); i++) {
 462                         if (m->control_pressed) { return 0; }
 463                         tool[i] = storage[i];
 464                 }
 465
 466                 for (int i = 0; i < row; i++){
 467                         if (m->control_pressed) { return 0; }
 468                         C1[i][0]=tool[i];
 469                         C1[i][1]=tool[i+row+row];
 470                         C1[i][2]=C1[i][1]/(secondGroupingStart);
 471
 472                         C2[i][0]=tool[i+row];
 473                         C2[i][1]=tool[i+row+row+row]; // var group 2
 474                         C2[i][2]=C2[i][1]/(column-secondGroupingStart);
 475                 }
 476
 477                 for (int i = 0; i < row; i++){
 478                         if (m->control_pressed) { return 0; }
 479                         xbardiff = C1[i][0]-C2[i][0];
 480                         denom = sqrt(C1[i][2]+C2[i][2]);
 481                         Ts[i]=fabs(xbardiff/denom);
 482                         if (fabs(Ts[i])>(fabs(Tinitial[i])+.0000000000001)){ //13th place
 483                                 counter[i]++;
 484                         }
 485                 }
 486
 487                 return 0;
 488
 489         }catch(exception& e) {
 490                 m->errorOut(e, "MothurMetastats", "calc_twosample_ts");
 491                 exit(1);
 492         }
 493 }
 494 /***********************************************************/
 495 vector<double> MothurMetastats::calc_qvalues(vector<double>& pValues) {
 496         try {
 497                 vector<double> qvalues;
 498                 int numRows = pValues.size();
 499
 500                 //fill lambdas with 0.00, 0.01, 0.02... 0.95
 501                 vector<double> lambdas(96, 0);
 502                 for (int i = 1; i < lambdas.size(); i++) { lambdas[i] = lambdas[i-1] + 0.01; }
 503
 504                 vector<double> pi0_hat(lambdas.size(), 0);
 505
 506                 //calculate pi0_hat
 507                 for (int l = 0; l < lambdas.size(); l++){ // for each lambda value
 508                         int count = 0;
 509                         for (int i = 0; i < numRows; i++){ // for each p-value in order
 510                                 if (pValues[i] > lambdas[l]){ count++; }
 511                                 pi0_hat[l] = count/(double)(numRows*(1-lambdas[l]));
 512                         }
 513                 }
 514
 515                 //from R code - replacing with spline and splint below
 516                 //vector<double> f_spline = smoothSpline(lambdas, pi0_hat, 3);
 517                 //double pi0 = f_spline[(f_spline.size()-1)];   // this is the essential pi0_hat value
 518
 519                 //cubic Spline
 520                 double pi0, notsure; //this is the essential pi0_hat value
 521                 int notSure;
 522                 vector<double> resultsSpline(lambdas.size(), 0.0);
 523                 spline(pi0_hat, lambdas, notSure, notSure, resultsSpline);
 524                 //some sort of loop to get best value??
 525                 splint(pi0_hat, lambdas, notsure, pi0, resultsSpline);
 526
 527                 //order p-values
 528                 vector<double> ordered_qs = qvalues;
 529                 vector<int> ordered_ps(pValues.size(), 0);
 530                 for (int i = 1; i < ordered_ps.size(); i++) {  ordered_ps[i] = ordered_ps[i-1] + 1; }
 531                 vector<double> tempPvalues = pValues;
 532                 OrderPValues(0, numRows-1, tempPvalues, ordered_ps);
 533
 534                 ordered_qs[numRows-1] <- min((pValues[ordered_ps[numRows-1]]*pi0), 1.0);
 535                 for (int i = (numRows-2); i >= 0; i--){
 536                         double p = pValues[ordered_ps[i]];
 537                         double temp = p*numRows*pi0/(double)i;
 538
 539                         ordered_qs[i] = min(temp, ordered_qs[i+1]);
 540                         ordered_qs[i] = min(ordered_qs[i], 1.0);
 541                 }
 542
 543                 //re-distribute calculated qvalues to appropriate rows
 544                 for (int i = 0; i < numRows; i++){
 545                         qvalues[ordered_ps[i]] = ordered_qs[i];
 546                 }
 547
 548                 return qvalues;
 549
 550         }catch(exception& e) {
 551                 m->errorOut(e, "MothurMetastats", "calc_qvalues");
 552                 exit(1);
 553         }
 554 }
 555 /***********************************************************/
 556 int MothurMetastats::OrderPValues(int low, int high, vector<double>& p, vector<int>& order) {
 557         try {
 558
 559                 if (low < high) {
 560                         int i = low+1;
 561                         int j = high;
 562                         int pivot = (low+high) / 2;
 563
 564                         swapElements(low, pivot, p, order);  //puts pivot in final spot
 565
 566                         /* compare value */
 567                         double key = p[low];
 568
 569                         /* partition */
 570                         while(i <= j) {
 571                                 /* find member above ... */
 572                                 while((i <= high) && (p[i] <= key))     {  i++;  }
 573
 574                                 /* find element below ... */
 575                                 while((j >= low) && (p[j] > key))       {  j--;  }
 576
 577                                 if(i < j) {
 578                                         swapElements(i, j, p, order);
 579                                 }
 580                         }
 581
 582                         swapElements(low, j, p, order);
 583
 584                         /* recurse */
 585                         OrderPValues(low, j-1, p, order);
 586                         OrderPValues(j+1, high, p, order);
 587                 }
 588
 589                 return 0;
 590
 591         }catch(exception& e) {
 592                 m->errorOut(e, "MothurMetastats", "OrderPValues");
 593                 exit(1);
 594         }
 595 }
 596 /***********************************************************/
 597 int MothurMetastats::swapElements(int i, int j, vector<double>& p, vector<int>& order) {
 598         try {
 599
 600                 double z = p[i];
 601                 p[i] = p[j];
 602                 p[j] = z;
 603
 604                 int temp = order[i];
 605                 order[i] = order[j];
 606                 order[j] = temp;
 607
 608                 return 0;
 609
 610         }catch(exception& e) {
 611                 m->errorOut(e, "MothurMetastats", "swapElements");
 612                 exit(1);
 613         }
 614 }
 615 /***********************************************************
 616 vector<double> MothurMetastats::smoothSpline(vector<double> x, vector<double> y, int df) {
 617         try {
 618
 619                 cout << "lambdas" << endl;
 620                 for (int l = 0; l < x.size(); l++){ cout << x[l] << '\t'; }
 621                 cout << endl << "pi0_hat" << endl;
 622                 for (int l = 0; l < y.size(); l++){ cout << y[l] << '\t'; }
 623                 cout << endl;
 624
 625                 //double low = -1.5; double high = 1.5; double tol = 1e-04; double eps = 2e-08; double maxit = 500;
 626                 int n = x.size();
 627                 vector<double> w(n, 1);
 628
 629                 //x <- signif(x, 6L) - I think this rounds to 6 decimals places
 630                 //ux <- unique(sort(x)) //x will be unique and sorted since we created it
 631                 //ox <- match(x, ux) //since its unique and sort it will match
 632
 633                 vector<double> wbar(n, 0);
 634                 vector<double> ybar(n, 0);
 635                 vector<double> xbar(n, 0.0);
 636                 double yssw = 0.0;
 637                 for (int i = 0; i < n; i++) {
 638                         wbar[i] = w[i];
 639                         ybar[i] = w[i]*y[i];
 640                         yssw += (w[i] * y[i] * y[i]) - wbar[i] * (ybar[i] * ybar[i]);
 641                         xbar[i] = (x[i] - x[0]) / (x[n-1] - x[0]);
 642                 }
 643
 644                 vector<double> knot = sknot1(xbar);
 645                 int nk = knot.size() - 4;
 646
 647                 //double ispar = 0.0; double spar = 0.0; double icrit = 3.0; double dofoff = 3.0;
 648
 649                 return y;
 650
 651         }catch(exception& e) {
 652                 m->errorOut(e, "MothurMetastats", "smoothSpline");
 653                 exit(1);
 654         }
 655 }
 656 /***********************************************************/
 657 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
 658 int MothurMetastats::spline(vector<double>& x, vector<double>& y, int yp1, int ypn, vector<double>& y2) {
 659         try {
 660
 661                 cout << "lambdas" << endl;
 662                 for (int l = 0; l < x.size(); l++){ cout << x[l] << '\t'; }
 663                 cout << endl << "pi0_hat" << endl;
 664                 for (int l = 0; l < y.size(); l++){ cout << y[l] << '\t'; }
 665                 cout << endl;
 666
 667                 double p, qn, sig, un;
 668
 669                 int n = y2.size();
 670                 vector<double> u(n-1, 0.0);
 671
 672                 if (yp1 > 0.99e30) { y2[0] = u[0] = 0.0; }
 673                 else {
 674                         y2[0] = -0.5;
 675                         u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1]-x[0]) - yp1);
 676                 }
 677
 678                 for (int i = 1; i < n-1; i++) {
 679                         sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]);
 680                         p = sig * y2[i-1] + 2.0;
 681                         y2[i] = (sig - 1.0)/p;
 682                         u[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i] - y[i-1]) / (x[i]-x[i-1]);
 683                         u[i] = (6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
 684                 }
 685
 686                 if (ypn > 0.99e30) { qn=un=0.0; }
 687                 else {
 688                         qn=0.5;
 689                         un=(3.0/(x[n-1]-x[n-2]))*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
 690                 }
 691
 692                 y2[n-1]=(un-qn*u[n-2])/(qn*y2[n-2]+1.0);
 693
 694                 for (int k=n-2; k>=0; k--) {
 695                         y2[k]=y2[k]*y2[k+1]+u[k];
 696                 }
 697
 698                 return 0;
 699
 700         }catch(exception& e) {
 701                 m->errorOut(e, "MothurMetastats", "spline");
 702                 exit(1);
 703         }
 704 }
 705 /***********************************************************/
 706 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
 707 int MothurMetastats::splint(vector<double>& xa, vector<double>& ya, double& x, double& y, vector<double>& y2a) {
 708         try {
 709                 int k;
 710                 double h,b,a;
 711
 712                 int n = xa.size();
 713
 714                 int klo=0;
 715                 int khi=n-1;
 716                 while (khi-klo > 1) {
 717
 718                         if (m->control_pressed) { break; }
 719
 720                         k = (khi+klo) >> 1;
 721                         if (xa[k] > x) { khi=k; }
 722                         else { klo=k; }
 723                 }
 724
 725                 h=xa[khi]-xa[klo];
 726                 if (h == 0.0) { m->mothurOut("[ERROR]: Bad xa input to splint routine."); m->mothurOutEndLine(); m->control_pressed = true; return 0; }
 727                 a=(xa[khi]-x)/h;
 728                 b=(x - xa[klo])/h;
 729                 y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
 730
 731                 return 0;
 732
 733         }catch(exception& e) {
 734                 m->errorOut(e, "MothurMetastats", "splint");
 735                 exit(1);
 736         }
 737 }
 738 /***********************************************************/
 739 vector<double> MothurMetastats::sknot1(vector<double> x) {
 740         try {
 741                 vector<double> knots;
 742                 int n = x.size();
 743                 int nk = nkn(n);
 744
 745                 cout << nk << endl;
 746                 //R equivalent - rep(x[1L], 3L)
 747                 knots.push_back(x[0]); knots.push_back(x[0]); knots.push_back(x[0]);
 748
 749                 //generate a sequence of nk equally spaced values from 1 to n. R equivalent = seq.int(1, n, length.out = nk)
 750                 vector<int> indexes = getSequence(0, n-1, nk);
 751                 for (int i = 0; i < indexes.size(); i++) { knots.push_back(x[indexes[i]]); }
 752
 753                 //R equivalent - rep(x[n], 3L)
 754                 knots.push_back(x[n-1]); knots.push_back(x[n-1]); knots.push_back(x[n-1]);
 755
 756                 return knots;
 757
 758         }catch(exception& e) {
 759                 m->errorOut(e, "MothurMetastats", "sknot1");
 760                 exit(1);
 761         }
 762 }
 763 /***********************************************************/
 764 vector<int> MothurMetastats::getSequence(int start, int end, int length) {
 765         try {
 766                 vector<int> sequence;
 767                 double increment = (end-start) / (double) (length-1);
 768
 769                 sequence.push_back(start);
 770                 for (int i = 1; i < length-1; i++) {
 771                         sequence.push_back(int(i*increment));
 772                 }
 773                 sequence.push_back(end);
 774
 775                 return sequence;
 776
 777         }catch(exception& e) {
 778                 m->errorOut(e, "MothurMetastats", "getSequence");
 779                 exit(1);
 780         }
 781 }
 782 /***********************************************************/
 783 //not right, havent fully figured out the variable types yet...
 784 int MothurMetastats::nkn(int n) {
 785         try {
 786
 787                 if (n < 50) { return n; }
 788                 else {
 789                         double a1 = log2(50);
 790                         double a2 = log2(100);
 791                         double a3 = log2(140);
 792                         double a4 = log2(200);
 793
 794                         if (n < 200) {
 795                                 return (int)pow(2.0, (a1 + (a2 - a1) * (n - 50)/(float)150));
 796                         }else if (n < 800) {
 797                                 return (int)pow(2.0, (a2 + (a3 - a2) * (n - 200)/(float)600));
 798                         }else if (n < 3200) {
 799                                 return (int)pow(2.0, (a3 + (a4 - a3) * (n - 800)/(float)2400));
 800                         }else {
 801                                 return (int)pow((double)(200 + (n - 3200)), 0.2);
 802                         }
 803                 }
 804
 805                 return 0;
 806
 807         }catch(exception& e) {
 808                 m->errorOut(e, "MothurMetastats", "nkn");
 809                 exit(1);
 810         }
 811 }
 812 /***********************************************************/
 813
 814
 815
 816
 817