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