trialSwap2.cpp

   1 #include "trialswap2.h"
   2
   3
   4 //The sum_of_squares, havel_hakimi and calc_c_score algorithms have been adapted from I. Miklos and J. Podani. 2004. Randomization of presence-absence matrices: comments and new algorithms. Ecology 85:86-92.
   5
   6
   7 /**************************************************************************************************/
   8 int TrialSwap2::intrand(int n){
   9     try {
  10         double z;
  11
  12         z = (double)random() * (double)n / (double)RAND_MAX;
  13         if(z>=n)
  14             z=n-1;
  15         if(z<0)
  16             z=0;
  17         return((int)floor(z));
  18     }
  19         catch(exception& e) {
  20                 m->errorOut(e, "TrialSwap2", "intrand");
  21                 exit(1);
  22         }
  23 }
  24 /**************************************************************************************************/
  25 /* completely random matrix, all column and row totals are variable, matrix size is the same
  26  *
  27  *
  28  */
  29 /**************************************************************************************************/
  30 int TrialSwap2::sim1(vector<vector<int> > &co_matrix){
  31     try {
  32         vector<int> randRow;
  33         vector<vector<int> > tmpmatrix;
  34         int nrows = co_matrix.size();
  35         int ncols = co_matrix[0].size();
  36
  37         //clear co_matrix
  38         //     for(i=0;i<nrows;i++)
  39         //     {
  40         //         co_matrix.clear();
  41         //     }
  42
  43         //cout << "building matrix" << endl;
  44         for(int i=0;i<nrows;i++){
  45             if (m->control_pressed) { break; }
  46
  47             for(int j=0;j<ncols;j++){
  48                 double randNum = rand() / double(RAND_MAX);
  49                 //cout << randNum << endl;
  50
  51                 if(randNum > 0.5) {
  52                     randRow.push_back(1);
  53                 }else{
  54                     randRow.push_back(0);
  55                 }
  56             }
  57             tmpmatrix.push_back(randRow);
  58             randRow.clear();
  59             //cout << endl;
  60         }
  61         co_matrix = tmpmatrix;
  62
  63         return 0;
  64     }
  65         catch(exception& e) {
  66                 m->errorOut(e, "TrialSwap2", "sim1");
  67                 exit(1);
  68         }
  69 }
  70 /**************************************************************************************************/
  71 /*
  72  *row sums fixed, columns equiprobable
  73  */
  74 void TrialSwap2::sim2(vector<vector<int> > &co_matrix)
  75 {
  76     try {
  77
  78         for(int i=0;i<co_matrix.size();i++)
  79         {
  80             if (m->control_pressed) { break; }
  81             random_shuffle( co_matrix[i].begin(), co_matrix[i].end() );
  82         }
  83     }
  84         catch(exception& e) {
  85                 m->errorOut(e, "TrialSwap2", "sim2");
  86                 exit(1);
  87         }
  88 }
  89 /**************************************************************************************************/
  90 int TrialSwap2::sim2plus(vector<int> rowtotal, vector<vector<int> > &co_matrix)
  91 {
  92     try {
  93         int nrows = co_matrix.size();
  94         int ncols = co_matrix[0].size();
  95         double cellprob = 1.0/ncols;
  96         vector<double> cellprobvec;
  97         vector<int> tmprow;
  98         vector<vector<int> > tmpmatrix;
  99         //double randNum;
 100
 101         double start = 0.0;
 102
 103         for(int i=0; i<ncols; i++)
 104         {
 105             if (m->control_pressed) { return 0; }
 106             cellprobvec.push_back(start + cellprob);
 107             start = cellprobvec[i];
 108         }
 109
 110         for(int i=0; i<nrows; i++)
 111         {
 112             tmprow.assign(ncols, 0);
 113
 114             while( accumulate( tmprow.begin(), tmprow.end(), 0 ) < rowtotal[i])
 115             {
 116                 if (m->control_pressed) { return 0; }
 117                 double randNum = rand() / double(RAND_MAX);
 118                 //cout << randNum << endl;
 119                 if(randNum <= cellprobvec[0])
 120                 {
 121                     tmprow[0] = 1;
 122                     continue;
 123                 }
 124                 for(int j=1;j<ncols;j++)
 125                 {
 126                     //cout << range[j] << endl;
 127                     if(randNum <= cellprobvec[j] && randNum > cellprobvec[j-1] && tmprow[j] != 1)
 128                     {
 129                         tmprow[j] = 1;
 130                     }
 131                 }
 132             }
 133             tmpmatrix.push_back(tmprow);
 134             tmprow.clear();
 135         }
 136         co_matrix = tmpmatrix;
 137         tmpmatrix.clear();
 138         cellprobvec.clear();
 139
 140         return 0;
 141     }
 142         catch(exception& e) {
 143                 m->errorOut(e, "TrialSwap2", "sim2plus");
 144                 exit(1);
 145         }
 146 }
 147 /**************************************************************************************************/
 148 /*
 149  * same as sim2 but using initmatrix which is the initial co-occurrence matrix before transposition
 150  * may have to be changed depending on what matrix 'seed' is used. One way to use is to transpose
 151  * every null matrix before using an index and use the random matrix as a seed for the next null.
 152  */
 153 /**************************************************************************************************/
 154 void TrialSwap2::sim3(vector<vector<int> > &initmatrix)
 155 {
 156     try {
 157         for(int i=0;i<initmatrix.size();i++)
 158         {
 159             if (m->control_pressed) { break; }
 160             random_shuffle( initmatrix[i].begin(), initmatrix[i].end() );
 161         }
 162
 163     }
 164         catch(exception& e) {
 165                 m->errorOut(e, "TrialSwap2", "sim3");
 166                 exit(1);
 167         }
 168 }
 169 /**************************************************************************************************/
 170 /*
 171  *
 172  *
 173  *
 174  */
 175 /**************************************************************************************************/
 176 int TrialSwap2::sim4(vector<int> columntotal, vector<int> rowtotal, vector<vector<int> > &co_matrix)
 177 {
 178     try {
 179         vector<double> colProb;
 180         vector<int> tmprow;//(ncols, 7);
 181         vector<vector<int> > tmpmatrix;
 182         vector<double> range;
 183         vector<double> randNums;
 184         int ncols = columntotal.size();
 185         int nrows = rowtotal.size();
 186         tmprow.clear();
 187
 188         double colSum = accumulate( columntotal.begin(), columntotal.end(), 0 );
 189         //cout << "col sum: " << colSum << endl;
 190         for(int i=0;i<ncols;i++)
 191         {
 192             if (m->control_pressed) { return 0; }
 193             colProb.push_back(columntotal[i]/colSum);
 194         }
 195
 196         double start = 0.0;
 197
 198         for(int i=0;i<ncols;i++)
 199         {
 200             if (m->control_pressed) { return 0; }
 201             range.push_back(start + colProb[i]);
 202             start = range[i];
 203         }
 204
 205         for(int i=0;i<nrows;i++)
 206         {
 207             tmprow.assign(ncols, 0);
 208             if (m->control_pressed) { return 0; }
 209
 210             while ( accumulate( tmprow.begin(), tmprow.end(), 0 ) < rowtotal[i])
 211             {
 212                 if (m->control_pressed) { return 0; }
 213
 214                 double randNum = rand() / double(RAND_MAX);
 215                 if(randNum <= range[0])
 216                 {
 217                     tmprow[0] = 1;
 218                     continue;
 219                 }
 220                 for(int j=1;j<ncols;j++)
 221                 {
 222                     if(randNum <= range[j] && randNum > range[j-1] && tmprow[j] != 1)
 223                     {
 224                         tmprow[j] = 1;
 225                     }
 226
 227                 }
 228             }
 229             tmpmatrix.push_back(tmprow);
 230             tmprow.clear();
 231         }
 232
 233         co_matrix = tmpmatrix;
 234
 235         return 0;
 236     }
 237         catch(exception& e) {
 238                 m->errorOut(e, "TrialSwap2", "sim4");
 239                 exit(1);
 240         }
 241 }
 242 /**************************************************************************************************/
 243 /*
 244  * inverse of sim4, MUST BE TRANSPOSED BEFORE CO-OCCURRENCE ANALYSIS
 245  *
 246  *
 247  */
 248 /**************************************************************************************************/
 249 int TrialSwap2::sim5(vector<int> initcolumntotal,vector<int> initrowtotal, vector<vector<int> > &initmatrix)
 250 {
 251     try {
 252         vector<double> colProb;
 253         vector<int> tmprow;//(ncols, 7);
 254         vector<vector<int> > tmpmatrix;
 255         vector<double> range;
 256         vector<double> randNums;
 257         int ncols = initcolumntotal.size();
 258         int nrows = initrowtotal.size();
 259
 260         tmprow.clear();
 261
 262         double colSum = accumulate( initcolumntotal.begin(), initcolumntotal.end(), 0 );
 263         //cout << "col sum: " << colSum << endl;
 264         for(int i=0;i<ncols;i++)
 265         {
 266             if (m->control_pressed) { return 0; }
 267             colProb.push_back(initcolumntotal[i]/colSum);
 268         }
 269
 270         double start = 0.0;
 271
 272         for(int i=0;i<ncols;i++)
 273         {
 274             if (m->control_pressed) { return 0; }
 275             range.push_back(start + colProb[i]);
 276             start = range[i];
 277         }
 278
 279         for(int i=0;i<nrows;i++)
 280         {
 281             tmprow.assign(ncols, 0);
 282             if (m->control_pressed) { return 0; }
 283
 284             while ( accumulate( tmprow.begin(), tmprow.end(), 0 ) < initrowtotal[i])
 285             {
 286                 if (m->control_pressed) { return 0; }
 287
 288                 double randNum = rand() / double(RAND_MAX);
 289                 if(randNum <= range[0])
 290                 {
 291                     tmprow[0] = 1;
 292                     continue;
 293                 }
 294                 for(int j=1;j<ncols;j++)
 295                 {
 296                     if(randNum <= range[j] && randNum > range[j-1] && tmprow[j] != 1)
 297                     {
 298                         tmprow[j] = 1;
 299                     }
 300
 301                 }
 302             }
 303             tmpmatrix.push_back(tmprow);
 304             tmprow.clear();
 305         }
 306
 307         initmatrix = tmpmatrix;
 308         return 0;
 309     }
 310         catch(exception& e) {
 311                 m->errorOut(e, "TrialSwap2", "sim5");
 312                 exit(1);
 313         }
 314 }
 315 /**************************************************************************************************/
 316 /*
 317  *
 318  *
 319  *
 320  */
 321 /**************************************************************************************************/
 322 int TrialSwap2::sim6(vector<int> columntotal, vector<vector<int> > &co_matrix)
 323 {
 324     try {
 325         vector<vector<int> > tmpmatrix;
 326         vector<double> colProb;
 327         vector<int> tmprow;
 328         vector<double> range;
 329         int ncols = columntotal.size();
 330         int nrows = co_matrix.size();
 331
 332         int colSum = accumulate( columntotal.begin(), columntotal.end(), 0 );
 333
 334         for(int i=0;i<ncols;i++)
 335         {
 336             if (m->control_pressed) { return 0; }
 337             colProb.push_back(columntotal[i]/double (colSum));
 338         }
 339
 340         double start = 0.0;
 341
 342         for(int i=0;i<ncols;i++)
 343         {
 344             if (m->control_pressed) { return 0; }
 345             range.push_back(start + colProb[i]);
 346             start = range[i];
 347         }
 348
 349         for(int i=0;i<nrows;i++)
 350         {
 351             if (m->control_pressed) { return 0; }
 352             tmprow.assign(ncols, 0);
 353             int tmprowtotal;
 354             tmprowtotal = (rand() / double (RAND_MAX)) * 10;
 355             while ( tmprowtotal > ncols) {
 356                 if (m->control_pressed) { return 0; }
 357                 tmprowtotal = (rand() / double (RAND_MAX)) * 10;
 358             }
 359             //cout << tmprowtotal << endl;
 360             //cout << accumulate( tmprow.begin(), tmprow.end(), 0 ) << endl;
 361
 362             while ( accumulate( tmprow.begin(), tmprow.end(), 0 ) < tmprowtotal)
 363             {
 364                 if (m->control_pressed) { return 0; }
 365                 double randNum = rand() / double(RAND_MAX);
 366                 //cout << randNum << endl;
 367                 if(randNum <= range[0])
 368                 {
 369                     tmprow[0] = 1;
 370                     continue;
 371                 }
 372                 for(int j=1;j<ncols;j++)
 373                 {
 374                     //cout << range[j] << endl;
 375                     if(randNum <= range[j] && randNum > range[j-1] && tmprow[j] != 1)
 376                     {
 377                         tmprow[j] = 1;
 378                     }
 379
 380                 }
 381
 382
 383             }
 384
 385             tmpmatrix.push_back(tmprow);
 386             tmprow.clear();
 387         }
 388
 389         co_matrix = tmpmatrix;
 390         tmpmatrix.clear();
 391
 392         return 0;
 393     }
 394         catch(exception& e) {
 395                 m->errorOut(e, "TrialSwap2", "sim6");
 396                 exit(1);
 397         }
 398 }
 399 /**************************************************************************************************/
 400 /*
 401  * MUST BE TRANSPOSED BEFORE CO-OCCURRENCE ANALYSIS
 402  *
 403  *
 404  */
 405 /**************************************************************************************************/
 406 int TrialSwap2::sim7(vector<int> initrowtotal, vector<vector<int> > &co_matrix)
 407 {
 408     try {
 409         vector<vector<double> > probmatrix;
 410         vector<vector<int> > tmpmatrix;
 411         vector<double> colProb;
 412         vector<double> probrow;
 413         vector<int> tmprow;
 414         vector<double> range;
 415         double nc;
 416         int ncols = co_matrix[0].size(); int nrows = co_matrix.size();
 417
 418         tmpmatrix.assign(nrows, vector<int>(ncols, 0.));
 419
 420         int rowsum = accumulate( initrowtotal.begin(), initrowtotal.end(), 0 );
 421
 422         nc = rowsum * ncols;
 423         //cout << nc << endl;
 424
 425         //assign null matrix based on probabilities
 426
 427         double start = 0.0; // don't reset start -- probs should be from 0-1 thoughout the entire matrix
 428
 429         for(int i=0;i<nrows;i++)
 430         {
 431             if (m->control_pressed) { return 0; }
 432             //cout << initrowtotal[i]/double(nc) << endl;
 433             double cellprob = initrowtotal[i]/double(nc);
 434             //cout << cellprob << endl;
 435             for(int j=0;j<ncols;j++)
 436             {
 437
 438                 probrow.push_back(start + cellprob);
 439                 //cout << probrow[j] << endl;
 440                 //cout << start << endl;
 441                 start = start + cellprob;
 442             }
 443             probmatrix.push_back(probrow);
 444             probrow.clear();
 445         }
 446
 447
 448         //while(tmprowsum < rowsum)
 449         //for(int k=0;k<rowsum;k++)
 450         int k = 0;
 451         while(k < rowsum)
 452         {
 453             if (m->control_pressed) { return 0; }
 454         done:
 455             //cout << k << endl;
 456             //tmprowsum = accumulate( tmprowtotal.begin(), tmprowtotal.end(), 0 );
 457             double randNum = rand() / double(RAND_MAX);
 458             //cout << randNum << "+" << endl;
 459             //special case for the first entry
 460             if(randNum <= probmatrix[0][0] && tmpmatrix[0][0] != 1)
 461             {
 462                 tmpmatrix[0][0] = 1;
 463                 k++;
 464                 //cout << k << endl;
 465                 continue;
 466             }
 467
 468
 469             for(int i=0;i<nrows;i++)
 470             {
 471                 if (m->control_pressed) { return 0; }
 472                 for(int j=0;j<ncols;j++)
 473                 {
 474                     //cout << probmatrix[i][j] << endl;
 475                     if(randNum <= probmatrix[i][j] && randNum > probmatrix[i][j-1] && tmpmatrix[i][j] != 1)
 476                     {
 477                         tmpmatrix[i][j] = 1;
 478                         k++;
 479                         //cout << k << endl;
 480                         goto done;
 481                     }
 482                     //else
 483                     //k = k-1;
 484                 }
 485
 486             }
 487
 488         }
 489
 490         co_matrix = tmpmatrix;
 491         return 0;
 492     //build probibility matrix
 493     /* for(int i=0;i<nrows;i++)
 494      {
 495      for(int j=0;j<ncols;j++)
 496      {
 497      probrow.push_back(rowtotal[i]/nc);
 498      }
 499      probmatrix.pushback(probrow);
 500      probrow.clear;
 501      }
 502      */
 503
 504     /* int colSum = accumulate( initcolumntotal.begin(), initcolumntotal.end(), 0 );
 505
 506         for(int i=0;i<ncols;i++)
 507         {
 508             colProb.push_back(initcolumntotal[i]/double (colSum));
 509         }
 510
 511         double start = 0.0;
 512
 513         for(int i=0;i<ncols;i++)
 514         {
 515             range.push_back(start + colProb[i]);
 516             start = range[i];
 517         }
 518
 519         for(int i=0;i<nrows;i++)
 520         {
 521             tmprow.assign(ncols, 0);
 522             int tmprowtotal;
 523             tmprowtotal = (rand() / double (RAND_MAX)) * 10;
 524             while ( tmprowtotal > ncols)
 525                 tmprowtotal = (rand() / double (RAND_MAX)) * 10;
 526             //cout << tmprowtotal << endl;
 527             //cout << accumulate( tmprow.begin(), tmprow.end(), 0 ) << endl;
 528
 529             while ( accumulate( tmprow.begin(), tmprow.end(), 0 ) < tmprowtotal)
 530             {
 531                 double randNum = rand() / double(RAND_MAX);
 532                 //cout << randNum << endl;
 533                 if(randNum <= range[0])
 534                 {
 535                     tmprow[0] = 1;
 536                     continue;
 537                 }
 538                 for(int j=1;j<ncols;j++)
 539                 {
 540                     //cout << range[j] << endl;
 541                     if(randNum <= range[j] && randNum > range[j-1] && tmprow[j] != 1)
 542                     {
 543                         tmprow[j] = 1;
 544                     }
 545                 }
 546             }
 547
 548             tmpmatrix.push_back(tmprow);
 549             tmprow.clear();
 550         }
 551
 552         initmatrix = tmpmatrix;
 553      */
 554     }
 555         catch(exception& e) {
 556                 m->errorOut(e, "TrialSwap2", "sim7");
 557                 exit(1);
 558         }
 559 }
 560 /**************************************************************************************************/
 561 /*
 562  *
 563  *
 564  *
 565  */
 566 /**************************************************************************************************/
 567 int TrialSwap2::sim8(vector<int> columntotal, vector<int> rowtotal, vector<vector<int> > &co_matrix)
 568 {
 569     try {
 570         double prob;
 571         double start = 0.0;
 572         int ncols = columntotal.size(); int nrows = rowtotal.size();
 573         double probarray[nrows * ncols];
 574         double randnum;
 575         int grandtotal;
 576         int total = 0;
 577
 578         //double colSum = accumulate( columntotal.begin(), columntotal.end(), 0 );
 579         double rowSum = accumulate( rowtotal.begin(), rowtotal.end(), 0 );
 580
 581         if (m->control_pressed) { return 0; }
 582
 583         //cout << "rowsum: " << rowSum << endl;
 584
 585         grandtotal = rowSum;
 586
 587         //create probability matrix with each site being between 0 and 1
 588         for (int i=0;i<nrows;i++) {
 589             if (m->control_pressed) { return 0; }
 590             for (int j=0;j<ncols;j++) {
 591                 prob = (rowtotal[i] * columntotal[j])/(rowSum*rowSum);
 592                 if (prob == 0.0)
 593                     probarray[ncols * i + j] = -1;
 594                 else
 595                     probarray[ncols * i + j] = start + prob;
 596                 //probmatrixrow.pushback(start + prob);
 597                 start += prob;
 598             }
 599         }
 600         //cout << "prbarray" << endl;
 601         //for(int i=0;i<(nrows*ncols);i++)
 602         //cout << probarray[i] << " ";
 603         //cout << endl;
 604
 605         //generate random muber between 0 and 1 and interate through probarray until found
 606         while (total < grandtotal)  {
 607             if (m->control_pressed) { return 0; }
 608             randnum = rand() / double(RAND_MAX);
 609             //cout << "rand num: " << randnum << endl;
 610             if((randnum <= probarray[0]) && (probarray[0] != 2) ) {
 611                 probarray[0] = 2;
 612                 total++;
 613                 continue;
 614             }
 615             for(int i=1;i<(nrows*ncols);i++) {
 616                 if (m->control_pressed) { return 0; }
 617                 if((randnum <= probarray[i]) && (randnum > probarray[i-1]) && (probarray[i] != 2) ) {
 618                     probarray[i] = 2;
 619                     total++;
 620                     break;
 621                 }
 622                 else
 623                     continue;
 624             }
 625         }
 626         //cout << "prbarray" << endl;
 627         //for(int i=0;i<(nrows*ncols);i++)
 628         //cout << probarray[i] << " ";
 629         //cout << endl;
 630         for(int i=0;i<nrows;i++) {
 631             if (m->control_pressed) { return 0; }
 632             for(int j=0;j<ncols;j++) {
 633                 if(probarray[ncols * i + j] == 2)
 634                     co_matrix[i][j] = 1;
 635                 else
 636                     co_matrix[i][j] = 0;
 637             }
 638         }
 639         return 0;
 640     }
 641         catch(exception& e) {
 642                 m->errorOut(e, "TrialSwap2", "sim8");
 643                 exit(1);
 644         }
 645 }
 646 /**************************************************************************************************/
 647 int TrialSwap2::havel_hakimi(vector<int> rowtotal,vector<int> columntotal,vector<vector<int> > &co_matrix)
 648 {
 649     try {
 650         int nrows = co_matrix.size();
 651         int ncols = co_matrix[0].size();
 652         int i,j,k;
 653         vector<int> r1; r1.resize(nrows,0);
 654         vector<int> c;  c.resize(ncols,0);
 655         vector<int> c1; c1.resize(ncols,0);
 656
 657
 658         for(i=0;i<nrows;i++) {
 659             for(j=0;j<ncols;j++) {
 660                 co_matrix[i][j]=0;
 661             }
 662         }
 663         for(i=0;i<nrows;i++) {
 664             r1[i]=1;
 665         }
 666
 667         for(i=0;i<ncols;i++)
 668         {
 669             c[i]=columntotal[i];
 670             c1[i]=i;
 671         }
 672
 673         for(k=0;k<nrows;k++)
 674         {
 675             if (m->control_pressed) { return 0; }
 676             i=intrand(nrows);
 677             while(r1[i]==0) {
 678                 if (m->control_pressed) { return 0; }
 679                 i=intrand(nrows);
 680             }
 681             r1[i]=0;
 682             sho(c,c1,ncols);
 683             for(j=0;j<rowtotal[i];j++)
 684             {
 685                 if (m->control_pressed) { return 0; }
 686                 co_matrix[i][c1[j]]=1;
 687                 c[j]--;
 688                 if(c[j]<0)
 689                     m->mothurOut("Uhh! " + toString(c1[j]) + "\n");
 690             }
 691         }
 692         return 0;
 693     }
 694         catch(exception& e) {
 695                 m->errorOut(e, "TrialSwap2", "havel_hakimi");
 696                 exit(1);
 697         }
 698 }
 699 /**************************************************************************************************/
 700 int TrialSwap2::sho(vector<int> c, vector<int> c1, int k)
 701 {
 702     try {
 703         int i,j,temp;
 704
 705         for(j=k-1;j>0;j--)
 706         {
 707             if (m->control_pressed) { return 0; }
 708             for(i=0;i<j;i++)
 709             {
 710                 if(c[i]<c[i+1])
 711                 {
 712                     temp=c[i];
 713                     c[i]=c[i+1];
 714                     c[i+1]=temp;
 715                     temp=c1[i];
 716                     c1[i]=c1[i+1];
 717                     c1[i+1]=temp;
 718                 }
 719             }
 720         }
 721         for(j=1;j<1000;j++)
 722         {
 723             if (m->control_pressed) { return 0; }
 724             i=intrand(k-1);
 725             if(c[i]==c[i+1])
 726             {
 727                 temp=c[i];
 728                 c[i]=c[i+1];
 729                 c[i+1]=temp;
 730                 temp=c1[i];
 731                 c1[i]=c1[i+1];
 732                 c1[i+1]=temp;
 733             }
 734         }
 735         return(0);
 736     }
 737         catch(exception& e) {
 738                 m->errorOut(e, "TrialSwap2", "sho");
 739                 exit(1);
 740         }
 741 }
 742 /**************************************************************************************************/
 743 double TrialSwap2::calc_c_score (vector<vector<int> > &co_matrix,vector<int>  rowtotal)
 744 {
 745     try {
 746         double cscore = 0.0;
 747         double maxD;
 748         double D;
 749         double normcscore = 0.0;
 750         int nonzeros = 0;
 751         int ncols = co_matrix[0].size(); int nrows = rowtotal.size();
 752         vector<vector<double> > s(nrows, vector<double>(nrows,0.0)); //only fill half the matrix
 753
 754         for(int i=0;i<nrows-1;i++)
 755         {
 756
 757             for(int j=i+1;j<nrows;j++)
 758             {
 759                 if (m->control_pressed) { return 0; }
 760                 for(int k=0;k<ncols;k++)
 761                 {
 762                     if((co_matrix[i][k]==1)&&(co_matrix[j][k]==1)) //if both are 1s ie co-occurrence
 763                         s[i][j]++; //s counts co-occurrences
 764                 }
 765
 766                 //rowtotal[i] = A, rowtotal[j] = B, ncols = P, s[i][j] = J
 767                 cscore += (rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);///(nrows*(nrows-1)/2);
 768                 D = (rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);
 769
 770                 if(ncols < (rowtotal[i] + rowtotal[j]))
 771                 {
 772                     maxD = (ncols-rowtotal[i])*(ncols-rowtotal[j]);
 773                 }
 774                 else
 775                 {
 776                     maxD = rowtotal[i] * rowtotal[j];
 777                 }
 778
 779                 if(maxD != 0)
 780                 {
 781                     normcscore += D/maxD;
 782                     nonzeros++;
 783                 }
 784             }
 785         }
 786
 787         cscore = cscore/(double)(nrows*(nrows-1)/2);
 788         //cout << "normalized c score: " << normcscore/nonzeros << endl;
 789
 790         return cscore;
 791     }
 792         catch(exception& e) {
 793                 m->errorOut(e, "TrialSwap2", "calc_c_score");
 794                 exit(1);
 795         }
 796 }
 797 /**************************************************************************************************/
 798 int TrialSwap2::calc_checker (vector<vector<int> > &co_matrix, vector<int>  rowtotal)
 799 {
 800     try {
 801         int cunits=0;
 802         //int s[nrows][ncols];
 803         int ncols = co_matrix[0].size(); int nrows = rowtotal.size();
 804         vector<vector<int> > s(nrows, vector<int>(nrows,0)); //only fill half the matrix
 805
 806         for(int i=0;i<nrows-1;i++)
 807         {
 808             for(int j=i+1;j<nrows;j++)
 809             {
 810                 if (m->control_pressed) { return 0; }
 811                 //s[i][j]=0;
 812                 for(int k=0;k<ncols;k++)
 813                 {
 814                     //cout << s[i][j] << endl;
 815                     //iterates through the row and counts co-occurrences. The total number of co-occurrences for each row pair is kept in matrix s at location s[i][j].
 816                     if((co_matrix[i][k]==1)&&(co_matrix[j][k]==1)) //if both are 1s ie co-occurrence
 817                         s[i][j]++; //s counts co-occurrences
 818
 819                 }
 820                 //cout << "rowtotal: " << rowtotal[i] << endl;
 821                 //cout << "co-occurrences: " << s[i][j] << endl;
 822                 //cunits+=(rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);
 823                 if (s[i][j] == 0)
 824                 {
 825                     cunits+=1;
 826                 }
 827                 //cunits+=s[i][j];
 828             }
 829         }
 830
 831         return cunits;
 832     }
 833         catch(exception& e) {
 834                 m->errorOut(e, "TrialSwap2", "calc_checker");
 835                 exit(1);
 836         }
 837 }
 838 /**************************************************************************************************/
 839 double TrialSwap2::calc_vratio (vector<int> rowtotal, vector<int> columntotal)
 840 {
 841     try {
 842         int nrows = rowtotal.size();
 843         int ncols = columntotal.size();
 844         int sumCol = accumulate(columntotal.begin(), columntotal.end(), 0 );
 845        // int sumRow = accumulate(rowtotal.begin(), rowtotal.end(), 0 );
 846
 847         double colAvg = (double) sumCol / (double) ncols;
 848  //       double rowAvg = (double) sumRow / (double) nrows;
 849
 850         double p = 0.0;
 851
 852  //       double totalRowVar = 0.0;
 853         double rowVar = 0.0;
 854         double colVar = 0.0;
 855
 856         for(int i=0;i<nrows;i++)
 857         {
 858             if (m->control_pressed) { return 0; }
 859             p = (double) rowtotal[i]/(double) ncols;
 860             rowVar += p * (1.0-p);
 861         }
 862
 863         for(int i=0;i<ncols;i++)
 864         {
 865             if (m->control_pressed) { return 0; }
 866             colVar += pow(((double) columntotal[i]-colAvg),2);
 867         }
 868
 869         colVar = (1.0/(double)ncols) * colVar;
 870
 871         return colVar/rowVar;
 872     }
 873     catch(exception& e) {
 874         m->errorOut(e, "TrialSwap2", "calc_vratio");
 875         exit(1);
 876     }
 877
 878 }
 879 /**************************************************************************************************/
 880 int TrialSwap2::calc_combo (vector<vector<int> > &initmatrix)
 881 {
 882     try {
 883         int initrows = initmatrix.size();
 884         int unique = 0;
 885         int match = 0;
 886         int matches = 0;
 887         for(int i=0;i<initrows;i++)
 888         {
 889             match = 0;
 890             for(int j=i+1;j<=initrows;j++)
 891             {
 892                 if (m->control_pressed) { return 0; }
 893                 if( (initmatrix[i] == initmatrix[j]))
 894                 {
 895                     match++;
 896                     matches++;
 897                     break;
 898                 }
 899             }
 900
 901             //on the last iteration of a previously matched row it will add itself because it doesn't match any following rows, so that combination is counted
 902             if (match == 0)
 903                 unique++;
 904         }
 905         return unique;
 906     }
 907     catch(exception& e) {
 908         m->errorOut(e, "TrialSwap2", "calc_combo");
 909         exit(1);
 910     }
 911 }
 912 /**************************************************************************************************/
 913 int TrialSwap2::swap_checkerboards (vector<vector<int> > &co_matrix)
 914 {
 915     try {
 916         int ncols = co_matrix[0].size(); int nrows = co_matrix.size();
 917         int i, j, k, l;
 918         i=intrand(nrows);
 919         while((j = intrand(nrows) ) == i ) {;if (m->control_pressed) { return 0; }}
 920         k=intrand(ncols);
 921         while((l = intrand(ncols) ) == k ) {;if (m->control_pressed) { return 0; }}
 922         //cout << co_matrix[i][k] << " " << co_matrix[j][l] << endl;
 923         //cout << co_matrix[i][l] << " " << co_matrix[j][k] << endl;
 924         //cout << co_matrix[i][l] << " " << co_matrix[j][k] << endl;
 925         //cout << co_matrix[i][l] << " " << co_matrix[j][k] << endl;
 926         if((co_matrix[i][k]*co_matrix[j][l]==1 && co_matrix[i][l]+co_matrix[j][k]==0)||(co_matrix[i][k]+co_matrix[j][l]==0 && co_matrix[i][l]*co_matrix[j][k]==1)) //checking for checkerboard value and swap
 927         {
 928             co_matrix[i][k]=1-co_matrix[i][k];
 929             co_matrix[i][l]=1-co_matrix[i][l];
 930             co_matrix[j][k]=1-co_matrix[j][k];
 931             co_matrix[j][l]=1-co_matrix[j][l];
 932             //cout << "swapped!" << endl;
 933         }
 934         //cout << "i: " << i << " j: " << j << " k: " << " l: " << l << endl;
 935         return 0;
 936     }
 937     catch(exception& e) {
 938         m->errorOut(e, "TrialSwap2", "swap_checkerboards");
 939         exit(1);
 940     }
 941 }
 942 /**************************************************************************************************/
 943 double TrialSwap2::calc_pvalue_greaterthan (vector<double> scorevec, double initialscore)
 944 {
 945     try {
 946         int runs = scorevec.size();
 947         double p = 0.0;
 948         for( int i=0;i<runs;i++)
 949         {
 950             if (m->control_pressed) { return 0; }
 951             if(scorevec[i]>=initialscore)
 952                 p++;
 953         }
 954         return p/(double)runs;
 955     }
 956     catch(exception& e) {
 957         m->errorOut(e, "TrialSwap2", "calc_pvalue_greaterthan");
 958         exit(1);
 959     }
 960 }
 961 /**************************************************************************************************/
 962 double TrialSwap2::calc_pvalue_lessthan (vector<double> scorevec, double initialscore)
 963 {
 964     try {
 965         int runs = scorevec.size();
 966         double p = 0.0;
 967         for( int i=0;i<runs;i++)
 968         {
 969             if (m->control_pressed) { return 0; }
 970             if(scorevec[i]<=initialscore)
 971                 p++;
 972         }
 973         return p/(double)runs;
 974     }
 975     catch(exception& e) {
 976         m->errorOut(e, "TrialSwap2", "calc_pvalue_lessthan");
 977         exit(1);
 978     }
 979 }
 980 /**************************************************************************************************/
 981 double TrialSwap2::t_test (double initialscore, int runs, double nullMean, vector<double> scorevec)
 982 {
 983     try {
 984         double t;
 985         double sampleSD;
 986         double sum = 0;
 987
 988         for(int i=0;i<runs;i++)
 989         {
 990             if (m->control_pressed) { return 0; }
 991             sum += pow((scorevec[i] - nullMean),2);
 992             //cout << "scorevec[" << i << "]" << scorevec[i] << endl;
 993         }
 994
 995         m->mothurOut("nullMean: " + toString(nullMean)); m->mothurOutEndLine();
 996
 997         m->mothurOut("sum: " + toString(sum));  m->mothurOutEndLine();
 998
 999         sampleSD = sqrt( (1/runs) * sum );
1000
1001         m->mothurOut("samplSD: " + toString(sampleSD));  m->mothurOutEndLine();
1002
1003         t = (nullMean - initialscore) / (sampleSD / sqrt(runs));
1004
1005         return t;
1006     }
1007     catch(exception& e) {
1008         m->errorOut(e, "TrialSwap2", "t_test");
1009         exit(1);
1010     }
1011 }
1012 /**************************************************************************************************/
1013 int TrialSwap2::print_matrix(vector<vector<int> > &matrix, int nrows, int ncols)
1014 {
1015     try {
1016          m->mothurOut("matrix:");  m->mothurOutEndLine();
1017
1018         for (int i = 0; i < nrows; i++)
1019         {
1020             if (m->control_pressed) { return 0; }
1021             for (int j = 0; j < ncols; j++)
1022             {
1023                 m->mothurOut(toString(matrix[i][j]));
1024             }
1025             m->mothurOutEndLine();
1026         }
1027         return 0;
1028     }
1029     catch(exception& e) {
1030         m->errorOut(e, "TrialSwap2", "print_matrix");
1031         exit(1);
1032     }
1033 }
1034 /**************************************************************************************************/
1035 int TrialSwap2::transpose_matrix (vector<vector<int> > &initmatrix, vector<vector<int> > &co_matrix)//, int nrows, int nocols)
1036 {
1037     try {
1038         int ncols = initmatrix.size(); int nrows = initmatrix[0].size();
1039         int tmpnrows = nrows;
1040         //vector<vector<int> > tmpvec;
1041         vector<int> tmprow;
1042         if(!co_matrix.empty())
1043             co_matrix.clear();
1044         for (int i=0;i<nrows;i++)
1045         {
1046             if (m->control_pressed) { return 0; }
1047             for (int j=0;j<ncols;j++)
1048             {
1049                 tmprow.push_back(initmatrix[j][i]);
1050             }
1051             /*if (accumulate( tmprow.begin(), tmprow.end(), 0 ) == 0)
1052              {
1053              tmpnrows--;
1054              }
1055              else */
1056             co_matrix.push_back(tmprow);
1057             tmprow.clear();
1058         }
1059         nrows = tmpnrows;
1060         return 0;
1061     }
1062     catch(exception& e) {
1063         m->errorOut(e, "TrialSwap2", "transpose_matrix");
1064         exit(1);
1065     }
1066 }
1067 /**************************************************************************************************/
1068 int TrialSwap2::update_row_col_totals(vector<vector<int> > &co_matrix, vector<int> &rowtotal, vector<int> &columntotal)
1069 {
1070     try {
1071         //rowtotal.clear();
1072         //columntotal.clear();
1073         //generate (rowtotal.begin(), rowtotal.end(), 0);
1074         //generate (columntotal.begin(), columntotal.end(), 0);
1075         int nrows = co_matrix.size();
1076         int ncols = co_matrix[0].size();
1077         vector<int> tmpcolumntotal(ncols, 0);
1078         vector<int> tmprowtotal(nrows, 0);
1079
1080         int rowcount = 0;
1081
1082         for (int i = 0; i < nrows; i++)
1083         {
1084             if (m->control_pressed) { return 0; }
1085             for (int j = 0; j < ncols; j++)
1086             {
1087                 if (co_matrix[i][j] == 1)
1088                 {
1089                     rowcount++;
1090                     tmpcolumntotal[j]++;
1091                 }
1092             }
1093             tmprowtotal[i] = rowcount;
1094             rowcount = 0;
1095         }
1096         columntotal = tmpcolumntotal;
1097         rowtotal = tmprowtotal;
1098         /*cout << "rowtotal: ";
1099         for(int i = 0; i<nrows; i++) { cout << rowtotal[i]; }
1100         cout << "  ";
1101         cout << " coltotal: ";
1102         for(int i = 0; i<ncols; i++) { cout << columntotal[i]; }
1103         cout << endl;*/
1104         return 0;
1105     }
1106     catch(exception& e) {
1107         m->errorOut(e, "TrialSwap2", "update_row_col_totals");
1108         exit(1);
1109     }
1110 }
1111 /**************************************************************************************************/
1112 /*int main(int argc, char *argv[])
1113 {
1114     srand (time(0));
1115     char* input_filename = argv[1];
1116     std::ifstream infile (input_filename);
1117
1118     //skip the first line of headers
1119     getline(infile, line);
1120     //get the first line of data
1121     getline(infile, line);
1122
1123     nrows = 0;
1124     ncols = 0;
1125
1126     //int numspaces = 0;
1127     char nextChar;
1128
1129     for (int i=0; i<int(line.length()); i++)
1130     {
1131       nextChar = line.at(i); // gets a character
1132       if (isspace(line[i]))
1133           ncols++;
1134     }
1135
1136     ncols = ncols-3;
1137
1138     cout << "number of OTUs: ";
1139     cout << ncols << endl;
1140
1141     infile.close();
1142
1143     std::ifstream infile2 (input_filename);
1144
1145     //skip first line of headers
1146     getline(infile2, line);
1147
1148     while (!infile2.eof())
1149     {
1150         getline(infile2, line);
1151         if (!line.empty())
1152             nrows++;
1153     }
1154
1155     cout << "number of sites: ";
1156     cout << nrows << endl;
1157
1158     infile2.close();
1159
1160     std::ifstream infile3 (input_filename);
1161
1162     //skip first line
1163     getline(infile3, line);
1164
1165     //variables that depend on info from initial matrix
1166     vector<vector<int> > co_matrix;//[nrows][ncols];
1167     vector<vector<int> > initmatrix;
1168     vector<int> tmprow;
1169     vector<double> stats;
1170     int tmpnrows = nrows;
1171
1172     for (int row1=0; row1<nrows; row1++) // first line was skipped when counting, so we can start from 0
1173     {
1174         //ignore first 3 cols in each row, data starts on the 3th col
1175         for (int i = 0; i < 3; i++)
1176             infile3 >> tmp;
1177
1178         for (int col=0; col<ncols; col++)
1179         {
1180             infile3 >> tmp;
1181             //cout << tmp << endl;
1182             if (atoi(tmp.c_str()) > 0)
1183                 tmprow.push_back(1);
1184             else
1185                 tmprow.push_back(0);
1186         }
1187         if (accumulate( tmprow.begin(), tmprow.end(), 0 ) == 0)
1188         {
1189             tmpnrows--;
1190         }
1191         else
1192             initmatrix.push_back(tmprow);
1193         //add the row to the matrix
1194         //initmatrix.push_back(tmprow);
1195         tmprow.clear();
1196         //cout << tmprow << endl;
1197     }
1198
1199     infile3.close();
1200     nrows = tmpnrows;
1201
1202     //print init matrix
1203     /* cout << "original matrix:" << endl;
1204
1205     for (int i = 0; i < nrows; i++)
1206     {
1207         for (int j = 0; j < ncols; j++)
1208         {
1209             cout << initmatrix[i][j];
1210         }
1211         cout << endl;
1212     } */
1213
1214         //for (i=0;i<ncols;i++)
1215         //cout << "col "<< i<< ": " << columntotal[i] << endl;
1216
1217     //co_matrix is now initmatrix and newmatrix is now co_matrix
1218
1219     //remove cols where sum is 0
1220
1221     //transpose matrix
1222    /* int newmatrows = ncols;
1223     int newmatcols = nrows;
1224     int initcols = ncols; //for the combo metric
1225     int initrows = nrows; //for the combo metric
1226     //swap for transposed matrix
1227     nrows = newmatrows;//ncols;
1228     ncols = newmatcols;//nrows;
1229
1230     vector<int> columntotal(ncols, 0);
1231     vector<int> initcolumntotal(ncols, 0);
1232     vector<int> initrowtotal(nrows, 0);
1233     vector<int> rowtotal(nrows, 0);
1234
1235     transpose_matrix(initmatrix,co_matrix);
1236     //remove degenerate rows and cols
1237
1238     //cout << "transposed matrix:" << endl;
1239     int rowcount = 0;
1240     for (int i = 0; i < nrows; i++)
1241     {
1242         for (int j = 0; j < ncols; j++)
1243         {
1244             if (co_matrix[i][j] == 1)
1245             {
1246                 rowcount++;
1247                 columntotal[j]++;
1248             }
1249             //cout << co_matrix[i][j];
1250         }
1251         //cout << " row total: " << rowcount << endl;
1252         //cout << endl;
1253         rowtotal[i] = rowcount;
1254         rowcount = 0;
1255     }
1256
1257     initcolumntotal = rowtotal;
1258     initrowtotal = columntotal;
1259
1260     cout << endl;
1261
1262     runs = atol(argv[2]);
1263     int metric = atol(argv[3]);
1264     int nullModel = atol(argv[4]);
1265     double initscore;
1266     update_row_col_totals(co_matrix, rowtotal, columntotal, ncols, nrows);
1267     //do initial metric: checker, c score, v ratio or combo
1268     switch(metric)
1269     {
1270         case 1:
1271             //c score
1272             initscore = calc_c_score(co_matrix, rowtotal);
1273             cout << "initial c score: " << initscore << endl;
1274             //print_matrix(co_matrix, nrows, ncols);
1275             break;
1276
1277         case 2:
1278             //checker
1279             initscore = calc_checker(co_matrix, rowtotal);
1280             cout << "initial checker score: " << initscore << endl;
1281             break;
1282
1283         case 3:
1284             //v ratio
1285             initscore = calc_vratio(nrows, ncols, rowtotal, columntotal);
1286             cout << "initial v ratio: " << initscore << endl;
1287             break;
1288
1289         case 4:
1290             //combo
1291             initscore = calc_combo(initrows, initcols, initmatrix);
1292             cout << "initial combo score: " << initscore << endl;
1293             //set co_matrix equal to initmatrix because combo requires row comparisons
1294             co_matrix = initmatrix;
1295             break;
1296
1297         case 5:
1298             //test!
1299
1300             //print_matrix(co_matrix, nrows, ncols);
1301             //sim1(nrows, ncols, co_matrix);
1302             //sim2(nrows, ncols, co_matrix);
1303             //sim3(initrows, initcols, initmatrix);
1304             //sim4(columntotal, rowtotal, co_matrix);
1305             //sim5(initcolumntotal, initmatrix);
1306             //sim6(columntotal, co_matrix);
1307             //sim7(initcolumntotal, initmatrix);
1308             sim8(columntotal, rowtotal, co_matrix);
1309             //print_matrix(initmatrix, initrows, initcols);
1310             //print_matrix(co_matrix, nrows, ncols);
1311
1312             break;
1313
1314         default:
1315             cout << "no metric selected!" << endl;
1316             return 1;
1317
1318     }
1319
1320     //matrix initialization
1321     //havel_hakimi(nrows, ncols, rowtotal, columntotal, co_matrix);
1322     //sum_of_square(nrows, ncols, rowtotal, columntotal, co_matrix);
1323     //co-matrix is now a random matrix with the same row and column totals as the initial matrix
1324
1325     //null matrix burn in
1326     cout << "initializing null matrix...";
1327     for(int l=0;l<10000;l++)
1328     {
1329        //swap_checkerboards (co_matrix);
1330        //if(l%10 == 0)
1331         switch(nullModel)
1332         {
1333             case 1:
1334                 //
1335                 sim1(nrows, ncols, co_matrix);
1336                 break;
1337
1338             case 2:
1339                 //sim2
1340                 sim2(nrows, ncols, co_matrix);
1341                 //sim2plus(nrows, ncols, initrowtotal, co_matrix);
1342                 break;
1343
1344             case 3:
1345                 //sim3
1346                 sim3(initrows, initcols, initmatrix);
1347                 //transpose_matrix(initmatrix,co_matrix);
1348                 co_matrix = initmatrix;
1349                 break;
1350
1351             case 4:
1352                 //sim4
1353                 sim4(columntotal, rowtotal, co_matrix);
1354                 break;
1355
1356             case 5:
1357                 //sim5
1358                 sim5(initcolumntotal, initrowtotal, initmatrix);
1359                 transpose_matrix(initmatrix,co_matrix);
1360                 //co_matrix = initmatrix;
1361                 break;
1362
1363             case 6:
1364                 sim6(columntotal, co_matrix);
1365                 break;
1366
1367             case 7:
1368                 //sim7(ncols, nrows, initrowtotal, co_matrix);
1369                 //transpose_matrix(initmatrix,co_matrix);
1370                 //co_matrix = initmatrix;
1371                 break;
1372
1373             case 8:
1374                 sim8(columntotal, rowtotal, co_matrix);
1375                 break;
1376
1377             case 9:
1378                 //swap_checkerboards
1379                 swap_checkerboards (co_matrix);
1380                 break;
1381
1382             default:
1383                 cout << "no null model selected!" << endl;
1384                 return 1;
1385         }
1386     }
1387     cout << "done!" << endl;
1388
1389     //generate null matrices and calculate the metrics
1390
1391     cout << "run: " << endl;
1392     for(int trial=0;trial<runs;trial++) //runs
1393     {
1394         printf("\b\b\b\b\b\b\b%7d",trial+1);
1395         fflush(stdout);
1396
1397         switch(nullModel)
1398         {
1399             case 1:
1400                 //
1401                 sim1(nrows, ncols, co_matrix);
1402                 break;
1403
1404             case 2:
1405                 //sim2
1406                 sim2(nrows, ncols, co_matrix);
1407                 //for(int i=0;i<nrows;i++)
1408                     //cout << rowtotal[i] << " ";
1409                 //sim2plus(nrows, ncols, initrowtotal, co_matrix);
1410                 break;
1411
1412             case 3:
1413                 //sim3
1414                 for(int i=0;i<nrows;i++)
1415                     cout  << " " << rowtotal[i];
1416                 sim3(initrows, initcols, initmatrix);
1417                 transpose_matrix(initmatrix,co_matrix);
1418                 break;
1419
1420             case 4:
1421                 //sim4
1422                 sim4(columntotal, rowtotal, co_matrix);
1423                 break;
1424
1425             case 5:
1426                 //sim5
1427                 sim5(initcolumntotal, initrowtotal, initmatrix);
1428                 transpose_matrix(initmatrix,co_matrix);
1429                 break;
1430
1431             case 6:
1432                 sim6(columntotal, co_matrix);
1433                 break;
1434
1435             case 7:
1436                 sim7(ncols, nrows, initrowtotal, co_matrix);
1437                 //print_matrix(co_matrix, nrows, ncols);
1438                 //transpose_matrix(initmatrix,co_matrix);
1439                 break;
1440
1441             case 8:
1442                 //sim8(initcolumntotal, initrowtotal, co_matrix);
1443                 //initrow and initcol are flipped because of transposition. this is super ugly!
1444                 sim8(initrowtotal, initcolumntotal, co_matrix);
1445                 break;
1446
1447             case 9:
1448                 //swap_checkerboards
1449                 swap_checkerboards (co_matrix);
1450                 break;
1451
1452             default:
1453                 cout << "no null model selected!" << endl;
1454                 return 1;
1455         }
1456
1457         //cout << endl;
1458         //print_matrix(co_matrix, nrows, ncols);
1459         update_row_col_totals(co_matrix, rowtotal, columntotal, ncols, nrows);
1460         //cout << columntotal[1]<<endl;
1461         double tmp;
1462         switch(metric)
1463         {
1464             case 1:
1465                 //c score
1466                 //swap_checkerboards(co_matrix);
1467                 //cout <<  "%" << tmp << " nrows " << nrows << " ncols " << ncols << " rowtotals ";
1468                 //for(int i = 0; i<nrows; i++) { cout << rowtotal[i]; }
1469                 //cout << endl;
1470                 tmp = calc_c_score(co_matrix, rowtotal);
1471                 //cout << "%" << tmp << " ";
1472                 stats.push_back(tmp);
1473                 //cout << "%" << tmp << " ";
1474                 //print_matrix(co_matrix, nrows, ncols);
1475                 break;
1476
1477             case 2:
1478                 //checker
1479                 //swap_checkerboards(co_matrix);
1480                 //cout <<  "%" << tmp << " nrows " << nrows << " ncols " << ncols << " rowtotals ";
1481                 //for(int i = 0; i<nrows; i++) { cout << rowtotal[i]; }
1482                 //cout << endl;
1483                 tmp = calc_checker(co_matrix, rowtotal);
1484                 stats.push_back(tmp);
1485                 //cout << "%" << tmp << endl;
1486                 break;
1487
1488             case 3:
1489                 //v ratio
1490                 stats.push_back(calc_vratio(nrows, ncols, rowtotal, columntotal));
1491                 break;
1492
1493             case 4:
1494                 //combo
1495                 stats.push_back(calc_combo(nrows, ncols, co_matrix));
1496                 break;
1497
1498             case 5:
1499                 //test!
1500                 break;
1501
1502             default:
1503                 cout << "no metric selected!" << endl;
1504                 return 1;
1505
1506         }
1507
1508     //print_matrix(co_matrix, nrows, ncols);
1509     //print_matrix(initmatrix, initrows, initcols);
1510
1511     }
1512
1513     cout << endl;
1514
1515     double total = 0.0;
1516     for (int i=0; i<stats.size();i++)
1517     {
1518         total+=stats[i];
1519     }
1520
1521     double nullMean = double (total/stats.size());
1522
1523     cout << '\n' << "average metric score: " << nullMean << endl;
1524
1525     //cout << "t-test: " << t_test(initscore, runs, nullMean, stats) << endl;
1526
1527     if (metric == 1 || metric == 2)
1528         cout << "pvalue: " << calc_pvalue_greaterthan (stats, initscore) << endl;
1529     else
1530         cout << "pvalue: " << calc_pvalue_lessthan (stats, initscore) << endl;
1531
1532     //print_matrix(co_matrix, nrows, ncols);
1533
1534     return 0;
1535
1536 }*/
1537 /**************************************************************************************************/
1538
1539