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::sim1(vector<vector<int> > &co_matrix){
   9     try {
  10         vector<int> randRow;
  11         vector<vector<int> > tmpmatrix;
  12         int nrows = co_matrix.size();
  13         int ncols = co_matrix[0].size();
  14
  15         //clear co_matrix
  16         //     for(i=0;i<nrows;i++)
  17         //     {
  18         //         co_matrix.clear();
  19         //     }
  20
  21         //cout << "building matrix" << endl;
  22         for(int i=0;i<nrows;i++){
  23             if (m->control_pressed) { break; }
  24
  25             for(int j=0;j<ncols;j++){
  26                 double randNum = rand() / double(RAND_MAX);
  27                 //cout << randNum << endl;
  28
  29                 if(randNum > 0.5) {
  30                     randRow.push_back(1);
  31                 }else{
  32                     randRow.push_back(0);
  33                 }
  34             }
  35             tmpmatrix.push_back(randRow);
  36             randRow.clear();
  37             //cout << endl;
  38         }
  39         co_matrix = tmpmatrix;
  40
  41         return 0;
  42     }
  43         catch(exception& e) {
  44                 m->errorOut(e, "TrialSwap2", "sim1");
  45                 exit(1);
  46         }
  47 }
  48 /**************************************************************************************************/
  49 /*
  50  *row sums fixed, columns equiprobable
  51  */
  52 void TrialSwap2::sim2(vector<vector<int> > &co_matrix)
  53 {
  54     try {
  55
  56         for(int i=0;i<co_matrix.size();i++)
  57         {
  58             if (m->control_pressed) { break; }
  59             random_shuffle( co_matrix[i].begin(), co_matrix[i].end() );
  60         }
  61     }
  62         catch(exception& e) {
  63                 m->errorOut(e, "TrialSwap2", "sim2");
  64                 exit(1);
  65         }
  66 }
  67 /**************************************************************************************************/
  68 int TrialSwap2::sim2plus(vector<int> rowtotal, vector<vector<int> > &co_matrix)
  69 {
  70     try {
  71         int nrows = co_matrix.size();
  72         int ncols = co_matrix[0].size();
  73         double cellprob = 1.0/ncols;
  74         vector<double> cellprobvec;
  75         vector<int> tmprow;
  76         vector<vector<int> > tmpmatrix;
  77         //double randNum;
  78
  79         double start = 0.0;
  80
  81         for(int i=0; i<ncols; i++)
  82         {
  83             if (m->control_pressed) { return 0; }
  84             cellprobvec.push_back(start + cellprob);
  85             start = cellprobvec[i];
  86         }
  87
  88         for(int i=0; i<nrows; i++)
  89         {
  90             tmprow.assign(ncols, 0);
  91
  92             while( accumulate( tmprow.begin(), tmprow.end(), 0 ) < rowtotal[i])
  93             {
  94                 if (m->control_pressed) { return 0; }
  95                 double randNum = rand() / double(RAND_MAX);
  96                 //cout << randNum << endl;
  97                 if(randNum <= cellprobvec[0])
  98                 {
  99                     tmprow[0] = 1;
 100                     continue;
 101                 }
 102                 for(int j=1;j<ncols;j++)
 103                 {
 104                     //cout << range[j] << endl;
 105                     if(randNum <= cellprobvec[j] && randNum > cellprobvec[j-1] && tmprow[j] != 1)
 106                     {
 107                         tmprow[j] = 1;
 108                     }
 109                 }
 110             }
 111             tmpmatrix.push_back(tmprow);
 112             tmprow.clear();
 113         }
 114         co_matrix = tmpmatrix;
 115         tmpmatrix.clear();
 116         cellprobvec.clear();
 117
 118         return 0;
 119     }
 120         catch(exception& e) {
 121                 m->errorOut(e, "TrialSwap2", "sim2plus");
 122                 exit(1);
 123         }
 124 }
 125 /**************************************************************************************************/
 126 /*
 127  * same as sim2 but using initmatrix which is the initial co-occurrence matrix before transposition
 128  * may have to be changed depending on what matrix 'seed' is used. One way to use is to transpose
 129  * every null matrix before using an index and use the random matrix as a seed for the next null.
 130  */
 131 /**************************************************************************************************/
 132 void TrialSwap2::sim3(vector<vector<int> > &initmatrix)
 133 {
 134     try {
 135         for(int i=0;i<initmatrix.size();i++)
 136         {
 137             if (m->control_pressed) { break; }
 138             random_shuffle( initmatrix[i].begin(), initmatrix[i].end() );
 139         }
 140
 141     }
 142         catch(exception& e) {
 143                 m->errorOut(e, "TrialSwap2", "sim3");
 144                 exit(1);
 145         }
 146 }
 147 /**************************************************************************************************/
 148 /*
 149  *
 150  *
 151  *
 152  */
 153 /**************************************************************************************************/
 154 int TrialSwap2::sim4(vector<int> columntotal, vector<int> rowtotal, vector<vector<int> > &co_matrix)
 155 {
 156     try {
 157         vector<double> colProb;
 158         vector<int> tmprow;//(ncols, 7);
 159         vector<vector<int> > tmpmatrix;
 160         vector<double> range;
 161         vector<double> randNums;
 162         int ncols = columntotal.size();
 163         int nrows = rowtotal.size();
 164         tmprow.clear();
 165
 166         double colSum = accumulate( columntotal.begin(), columntotal.end(), 0 );
 167         //cout << "col sum: " << colSum << endl;
 168         for(int i=0;i<ncols;i++)
 169         {
 170             if (m->control_pressed) { return 0; }
 171             colProb.push_back(columntotal[i]/colSum);
 172         }
 173
 174         double start = 0.0;
 175
 176         for(int i=0;i<ncols;i++)
 177         {
 178             if (m->control_pressed) { return 0; }
 179             range.push_back(start + colProb[i]);
 180             start = range[i];
 181         }
 182
 183         for(int i=0;i<nrows;i++)
 184         {
 185             tmprow.assign(ncols, 0);
 186             if (m->control_pressed) { return 0; }
 187
 188             while ( accumulate( tmprow.begin(), tmprow.end(), 0 ) < rowtotal[i])
 189             {
 190                 if (m->control_pressed) { return 0; }
 191
 192                 double randNum = rand() / double(RAND_MAX);
 193                 if(randNum <= range[0])
 194                 {
 195                     tmprow[0] = 1;
 196                     continue;
 197                 }
 198                 for(int j=1;j<ncols;j++)
 199                 {
 200                     if(randNum <= range[j] && randNum > range[j-1] && tmprow[j] != 1)
 201                     {
 202                         tmprow[j] = 1;
 203                     }
 204
 205                 }
 206             }
 207             tmpmatrix.push_back(tmprow);
 208             tmprow.clear();
 209         }
 210
 211         co_matrix = tmpmatrix;
 212
 213         return 0;
 214     }
 215         catch(exception& e) {
 216                 m->errorOut(e, "TrialSwap2", "sim4");
 217                 exit(1);
 218         }
 219 }
 220 /**************************************************************************************************/
 221 /*
 222  * inverse of sim4, MUST BE TRANSPOSED BEFORE CO-OCCURRENCE ANALYSIS
 223  *
 224  *
 225  */
 226 /**************************************************************************************************/
 227 int TrialSwap2::sim5(vector<int> initcolumntotal,vector<int> initrowtotal, vector<vector<int> > &initmatrix)
 228 {
 229     try {
 230         vector<double> colProb;
 231         vector<int> tmprow;//(ncols, 7);
 232         vector<vector<int> > tmpmatrix;
 233         vector<double> range;
 234         vector<double> randNums;
 235         int ncols = initcolumntotal.size();
 236         int nrows = initrowtotal.size();
 237
 238         tmprow.clear();
 239
 240         double colSum = accumulate( initcolumntotal.begin(), initcolumntotal.end(), 0 );
 241         //cout << "col sum: " << colSum << endl;
 242         for(int i=0;i<ncols;i++)
 243         {
 244             if (m->control_pressed) { return 0; }
 245             colProb.push_back(initcolumntotal[i]/colSum);
 246         }
 247
 248         double start = 0.0;
 249
 250         for(int i=0;i<ncols;i++)
 251         {
 252             if (m->control_pressed) { return 0; }
 253             range.push_back(start + colProb[i]);
 254             start = range[i];
 255         }
 256
 257         for(int i=0;i<nrows;i++)
 258         {
 259             tmprow.assign(ncols, 0);
 260             if (m->control_pressed) { return 0; }
 261
 262             while ( accumulate( tmprow.begin(), tmprow.end(), 0 ) < initrowtotal[i])
 263             {
 264                 if (m->control_pressed) { return 0; }
 265
 266                 double randNum = rand() / double(RAND_MAX);
 267                 if(randNum <= range[0])
 268                 {
 269                     tmprow[0] = 1;
 270                     continue;
 271                 }
 272                 for(int j=1;j<ncols;j++)
 273                 {
 274                     if(randNum <= range[j] && randNum > range[j-1] && tmprow[j] != 1)
 275                     {
 276                         tmprow[j] = 1;
 277                     }
 278
 279                 }
 280             }
 281             tmpmatrix.push_back(tmprow);
 282             tmprow.clear();
 283         }
 284
 285         initmatrix = tmpmatrix;
 286         return 0;
 287     }
 288         catch(exception& e) {
 289                 m->errorOut(e, "TrialSwap2", "sim5");
 290                 exit(1);
 291         }
 292 }
 293 /**************************************************************************************************/
 294 /*
 295  *
 296  *
 297  *
 298  */
 299 /**************************************************************************************************/
 300 int TrialSwap2::sim6(vector<int> columntotal, vector<vector<int> > &co_matrix)
 301 {
 302     try {
 303         vector<vector<int> > tmpmatrix;
 304         vector<double> colProb;
 305         vector<int> tmprow;
 306         vector<double> range;
 307         int ncols = columntotal.size();
 308         int nrows = co_matrix.size();
 309
 310         int colSum = accumulate( columntotal.begin(), columntotal.end(), 0 );
 311
 312         for(int i=0;i<ncols;i++)
 313         {
 314             if (m->control_pressed) { return 0; }
 315             colProb.push_back(columntotal[i]/double (colSum));
 316         }
 317
 318         double start = 0.0;
 319
 320         for(int i=0;i<ncols;i++)
 321         {
 322             if (m->control_pressed) { return 0; }
 323             range.push_back(start + colProb[i]);
 324             start = range[i];
 325         }
 326
 327         for(int i=0;i<nrows;i++)
 328         {
 329             if (m->control_pressed) { return 0; }
 330             tmprow.assign(ncols, 0);
 331             int tmprowtotal;
 332             tmprowtotal = (rand() / double (RAND_MAX)) * 10;
 333             while ( tmprowtotal > ncols) {
 334                 if (m->control_pressed) { return 0; }
 335                 tmprowtotal = (rand() / double (RAND_MAX)) * 10;
 336             }
 337             //cout << tmprowtotal << endl;
 338             //cout << accumulate( tmprow.begin(), tmprow.end(), 0 ) << endl;
 339
 340             while ( accumulate( tmprow.begin(), tmprow.end(), 0 ) < tmprowtotal)
 341             {
 342                 if (m->control_pressed) { return 0; }
 343                 double randNum = rand() / double(RAND_MAX);
 344                 //cout << randNum << endl;
 345                 if(randNum <= range[0])
 346                 {
 347                     tmprow[0] = 1;
 348                     continue;
 349                 }
 350                 for(int j=1;j<ncols;j++)
 351                 {
 352                     //cout << range[j] << endl;
 353                     if(randNum <= range[j] && randNum > range[j-1] && tmprow[j] != 1)
 354                     {
 355                         tmprow[j] = 1;
 356                     }
 357
 358                 }
 359
 360
 361             }
 362
 363             tmpmatrix.push_back(tmprow);
 364             tmprow.clear();
 365         }
 366
 367         co_matrix = tmpmatrix;
 368         tmpmatrix.clear();
 369
 370         return 0;
 371     }
 372         catch(exception& e) {
 373                 m->errorOut(e, "TrialSwap2", "sim6");
 374                 exit(1);
 375         }
 376 }
 377 /**************************************************************************************************/
 378 /*
 379  * MUST BE TRANSPOSED BEFORE CO-OCCURRENCE ANALYSIS
 380  *
 381  *
 382  */
 383 /**************************************************************************************************/
 384 int TrialSwap2::sim7(vector<int> initrowtotal, vector<vector<int> > &co_matrix)
 385 {
 386     try {
 387         vector<vector<double> > probmatrix;
 388         vector<vector<int> > tmpmatrix;
 389         vector<double> colProb;
 390         vector<double> probrow;
 391         vector<int> tmprow;
 392         vector<double> range;
 393         double nc;
 394         int ncols = co_matrix[0].size(); int nrows = co_matrix.size();
 395
 396         tmpmatrix.assign(nrows, vector<int>(ncols, 0.));
 397
 398         int rowsum = accumulate( initrowtotal.begin(), initrowtotal.end(), 0 );
 399
 400         nc = rowsum * ncols;
 401         //cout << nc << endl;
 402
 403         //assign null matrix based on probabilities
 404
 405         double start = 0.0; // don't reset start -- probs should be from 0-1 thoughout the entire matrix
 406
 407         for(int i=0;i<nrows;i++)
 408         {
 409             if (m->control_pressed) { return 0; }
 410             //cout << initrowtotal[i]/double(nc) << endl;
 411             double cellprob = initrowtotal[i]/double(nc);
 412             //cout << cellprob << endl;
 413             for(int j=0;j<ncols;j++)
 414             {
 415
 416                 probrow.push_back(start + cellprob);
 417                 //cout << probrow[j] << endl;
 418                 //cout << start << endl;
 419                 start = start + cellprob;
 420             }
 421             probmatrix.push_back(probrow);
 422             probrow.clear();
 423         }
 424
 425
 426         //while(tmprowsum < rowsum)
 427         //for(int k=0;k<rowsum;k++)
 428         int k = 0;
 429         while(k < rowsum)
 430         {
 431             if (m->control_pressed) { return 0; }
 432         done:
 433             //cout << k << endl;
 434             //tmprowsum = accumulate( tmprowtotal.begin(), tmprowtotal.end(), 0 );
 435             double randNum = rand() / double(RAND_MAX);
 436             //cout << randNum << "+" << endl;
 437             //special case for the first entry
 438             if(randNum <= probmatrix[0][0] && tmpmatrix[0][0] != 1)
 439             {
 440                 tmpmatrix[0][0] = 1;
 441                 k++;
 442                 //cout << k << endl;
 443                 continue;
 444             }
 445
 446
 447             for(int i=0;i<nrows;i++)
 448             {
 449                 if (m->control_pressed) { return 0; }
 450                 for(int j=0;j<ncols;j++)
 451                 {
 452                     //cout << probmatrix[i][j] << endl;
 453                     if(randNum <= probmatrix[i][j] && randNum > probmatrix[i][j-1] && tmpmatrix[i][j] != 1)
 454                     {
 455                         tmpmatrix[i][j] = 1;
 456                         k++;
 457                         //cout << k << endl;
 458                         goto done;
 459                     }
 460                     //else
 461                     //k = k-1;
 462                 }
 463
 464             }
 465
 466         }
 467
 468         co_matrix = tmpmatrix;
 469         return 0;
 470     //build probibility matrix
 471     /* for(int i=0;i<nrows;i++)
 472      {
 473      for(int j=0;j<ncols;j++)
 474      {
 475      probrow.push_back(rowtotal[i]/nc);
 476      }
 477      probmatrix.pushback(probrow);
 478      probrow.clear;
 479      }
 480      */
 481
 482     /* int colSum = accumulate( initcolumntotal.begin(), initcolumntotal.end(), 0 );
 483
 484         for(int i=0;i<ncols;i++)
 485         {
 486             colProb.push_back(initcolumntotal[i]/double (colSum));
 487         }
 488
 489         double start = 0.0;
 490
 491         for(int i=0;i<ncols;i++)
 492         {
 493             range.push_back(start + colProb[i]);
 494             start = range[i];
 495         }
 496
 497         for(int i=0;i<nrows;i++)
 498         {
 499             tmprow.assign(ncols, 0);
 500             int tmprowtotal;
 501             tmprowtotal = (rand() / double (RAND_MAX)) * 10;
 502             while ( tmprowtotal > ncols)
 503                 tmprowtotal = (rand() / double (RAND_MAX)) * 10;
 504             //cout << tmprowtotal << endl;
 505             //cout << accumulate( tmprow.begin(), tmprow.end(), 0 ) << endl;
 506
 507             while ( accumulate( tmprow.begin(), tmprow.end(), 0 ) < tmprowtotal)
 508             {
 509                 double randNum = rand() / double(RAND_MAX);
 510                 //cout << randNum << endl;
 511                 if(randNum <= range[0])
 512                 {
 513                     tmprow[0] = 1;
 514                     continue;
 515                 }
 516                 for(int j=1;j<ncols;j++)
 517                 {
 518                     //cout << range[j] << endl;
 519                     if(randNum <= range[j] && randNum > range[j-1] && tmprow[j] != 1)
 520                     {
 521                         tmprow[j] = 1;
 522                     }
 523                 }
 524             }
 525
 526             tmpmatrix.push_back(tmprow);
 527             tmprow.clear();
 528         }
 529
 530 /**************************************************************************************************/
 531 double TrialSwap2::calc_c_score (vector<vector<int> > &co_matrix, vector<int>  rowtotal, int ncols, int nrows)
 532 {
 533     try {
 534         double cscore = 0.0;
 535         double maxD;
 536         double D;
 537         double normcscore = 0.0;
 538         int nonzeros = 0;
 539         //int ncols = co_matrix[0].size(); int nrows = rowtotal.size();
 540         vector<vector<double> > s; s.resize(nrows);
 541         for (int i = 0; i < nrows; i++) { s[i].resize(nrows,0.0); }//only fill half the matrix
 542
 543
 544         for(int i=0;i<nrows-1;i++)
 545         {
 546
 547             for(int j=i+1;j<nrows;j++)
 548             {
 549                 if (m->control_pressed) { return 0; }
 550                 for(int k=0;k<ncols;k++)
 551                 {
 552                     if((co_matrix[i][k]==1)&&(co_matrix[j][k]==1)) //if both are 1s ie co-occurrence
 553                         s[i][j]++; //s counts co-occurrences
 554                 }
 555
 556                 //rowtotal[i] = A, rowtotal[j] = B, ncols = P, s[i][j] = J
 557                 cscore += (rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);///(nrows*(nrows-1)/2);
 558                 D = (rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);
 559
 560                 if(ncols < (rowtotal[i] + rowtotal[j]))
 561                 {
 562                     maxD = (ncols-rowtotal[i])*(ncols-rowtotal[j]);
 563                 }
 564                 else
 565                 {
 566                     maxD = rowtotal[i] * rowtotal[j];
 567                 }
 568
 569                 if(maxD != 0)
 570                 {
 571                     normcscore += D/maxD;
 572                     nonzeros++;
 573                 }
 574             }
 575         }
 576
 577         cscore = cscore/(double)(nrows*(nrows-1)/2);
 578         //cout << "normalized c score: " << normcscore/nonzeros << endl;
 579
 580         return cscore;
 581     }
 582         catch(exception& e) {
 583                 m->errorOut(e, "TrialSwap2", "calc_c_score");
 584                 exit(1);
 585         }
 586 }
 587 /**************************************************************************************************/
 588 int TrialSwap2::calc_checker (vector<vector<int> > &co_matrix, vector<int>  rowtotal, int ncols, int nrows)
 589 {
 590     try {
 591         int cunits=0;
 592         //int s[nrows][ncols];
 593         //int ncols = co_matrix[0].size(); int nrows = rowtotal.size();
 594         vector<vector<int> > s; s.resize(nrows);
 595         for (int i = 0; i < nrows; i++) { s[i].resize(nrows,0); }//only fill half the matrix
 596
 597         for(int i=0;i<nrows-1;i++)
 598         {
 599             for(int j=i+1;j<nrows;j++)
 600             {
 601                 if (m->control_pressed) { return 0; }
 602                 //s[i][j]=0;
 603                 for(int k=0;k<ncols;k++)
 604                 {
 605                     //cout << s[i][j] << endl;
 606                     //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].
 607                     if((co_matrix[i][k]==1)&&(co_matrix[j][k]==1)) //if both are 1s ie co-occurrence
 608                         s[i][j]++; //s counts co-occurrences
 609
 610                 }
 611                 //cout << "rowtotal: " << rowtotal[i] << endl;
 612                 //cout << "co-occurrences: " << s[i][j] << endl;
 613                 //cunits+=(rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);
 614                 if (s[i][j] == 0)
 615                 {
 616                     cunits+=1;
 617                 }
 618                 //cunits+=s[i][j];
 619             }
 620         }
 621
 622         return cunits;
 623     }
 624         catch(exception& e) {
 625                 m->errorOut(e, "TrialSwap2", "calc_checker");
 626                 exit(1);
 627         }
 628 }
 629 /**************************************************************************************************/
 630 double TrialSwap2::calc_vratio (int nrows, int ncols, vector<int> rowtotal, vector<int> columntotal)
 631 {
 632     try {
 633         //int nrows = rowtotal.size();
 634         //int ncols = columntotal.size();
 635         int sumCol = accumulate(columntotal.begin(), columntotal.end(), 0 );
 636        // int sumRow = accumulate(rowtotal.begin(), rowtotal.end(), 0 );
 637
 638         double colAvg = (double) sumCol / (double) ncols;
 639  //       double rowAvg = (double) sumRow / (double) nrows;
 640
 641         double p = 0.0;
 642
 643  //       double totalRowVar = 0.0;
 644         double rowVar = 0.0;
 645         double colVar = 0.0;
 646
 647         for(int i=0;i<nrows;i++)
 648         {
 649             if (m->control_pressed) { return 0; }
 650             p = (double) rowtotal[i]/(double) ncols;
 651             rowVar += p * (1.0-p);
 652         }
 653
 654         for(int i=0;i<ncols;i++)
 655         {
 656             if (m->control_pressed) { return 0; }
 657             colVar += pow(((double) columntotal[i]-colAvg),2);
 658         }
 659
 660         colVar = (1.0/(double)ncols) * colVar;
 661
 662         return colVar/rowVar;
 663     }
 664     catch(exception& e) {
 665         m->errorOut(e, "TrialSwap2", "calc_vratio");
 666         exit(1);
 667     }
 668
 669 }
 670 /**************************************************************************************************/
 671 int TrialSwap2::calc_combo (int nrows, int ncols, vector<vector<int> > &nullmatrix)
 672 {
 673     try {
 674         //need to transpose so we can compare rows (row-major order)
 675         int tmpnrows = nrows;
 676         vector<vector<int> > tmpmatrix;
 677
 678         vector<int> tmprow;
 679         if(!tmpmatrix.empty())
 680             tmpmatrix.clear();
 681         for (int i=0;i<ncols;i++)
 682         {
 683             for (int j=0;j<nrows;j++)
 684             {
 685                 tmprow.push_back(nullmatrix[j][i]);
 686             }
 687
 688             tmpmatrix.push_back(tmprow);
 689             tmprow.clear();
 690         }
 691
 692         int unique = 0;
 693         int match = 0;
 694         for(int j=0;j<ncols;j++)
 695         {
 696             match = 0;
 697             for(int i=j+1;i<=ncols;i++)
 698             {
 699                 //comparing matrix rows
 700                 if( (tmpmatrix[j] == tmpmatrix[i]))
 701                 {
 702                     match++;
 703                     break;
 704                 }
 705             }
 706
 707             //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
 708             if (match == 0)
 709                 unique++;
 710     }
 711     return unique;
 712     }
 713     catch(exception& e) {
 714         m->errorOut(e, "TrialSwap2", "calc_combo");
 715         exit(1);
 716     }
 717 }
 718 /**************************************************************************************************/
 719 int TrialSwap2::swap_checkerboards (vector<vector<int> > &co_matrix)
 720 {
 721     try {
 722         int ncols = co_matrix[0].size(); int nrows = co_matrix.size();
 723         int i, j, k, l;
 724         i = m->getRandomIndex(nrows-1);
 725         while((j = m->getRandomIndex(nrows-1) ) == i ) {;if (m->control_pressed) { return 0; }}
 726         k = m->getRandomIndex(ncols-1);
 727         while((l = m->getRandomIndex(ncols-1)) == k ) {;if (m->control_pressed) { return 0; }}
 728
 729         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
 730         {
 731             co_matrix[i][k]=1-co_matrix[i][k];
 732             co_matrix[i][l]=1-co_matrix[i][l];
 733             co_matrix[j][k]=1-co_matrix[j][k];
 734             co_matrix[j][l]=1-co_matrix[j][l];
 735
 736         }
 737
 738         return 0;
 739     }
 740     catch(exception& e) {
 741         m->errorOut(e, "TrialSwap2", "swap_checkerboards");
 742         exit(1);
 743     }
 744 }
 745 /**************************************************************************************************/
 746 double TrialSwap2::calc_pvalue_greaterthan (vector<double> scorevec, double initialscore)
 747 {
 748     try {
 749         int runs = scorevec.size();
 750         double p = 0.0;
 751         for( int i=0;i<runs;i++)
 752         {
 753             if (m->control_pressed) { return 0; }
 754             if(scorevec[i]>=initialscore)
 755                 p++;
 756         }
 757         return p/(double)runs;
 758     }
 759     catch(exception& e) {
 760         m->errorOut(e, "TrialSwap2", "calc_pvalue_greaterthan");
 761         exit(1);
 762     }
 763 }
 764 /**************************************************************************************************/
 765 double TrialSwap2::calc_pvalue_lessthan (vector<double> scorevec, double initialscore)
 766 {
 767     try {
 768         int runs = scorevec.size();
 769         double p = 0.0;
 770         for( int i=0;i<runs;i++)
 771         {
 772             if (m->control_pressed) { return 0; }
 773             if(scorevec[i]<=initialscore)
 774                 p++;
 775         }
 776         return p/(double)runs;
 777     }
 778     catch(exception& e) {
 779         m->errorOut(e, "TrialSwap2", "calc_pvalue_lessthan");
 780         exit(1);
 781     }
 782 }
 783 /**************************************************************************************************/
 784 double TrialSwap2::t_test (double initialscore, int runs, double nullMean, vector<double> scorevec)
 785 {
 786     try {
 787         double t;
 788         double sampleSD;
 789         double sum = 0;
 790
 791         for(int i=0;i<runs;i++)
 792         {
 793             if (m->control_pressed) { return 0; }
 794             sum += pow((scorevec[i] - nullMean),2);
 795             //cout << "scorevec[" << i << "]" << scorevec[i] << endl;
 796         }
 797
 798         m->mothurOut("nullMean: " + toString(nullMean)); m->mothurOutEndLine();
 799
 800         m->mothurOut("sum: " + toString(sum));  m->mothurOutEndLine();
 801
 802         sampleSD = sqrt( (1/runs) * sum );
 803
 804         m->mothurOut("samplSD: " + toString(sampleSD));  m->mothurOutEndLine();
 805
 806         t = (nullMean - initialscore) / (sampleSD / sqrt(runs));
 807
 808         return t;
 809     }
 810     catch(exception& e) {
 811         m->errorOut(e, "TrialSwap2", "t_test");
 812         exit(1);
 813     }
 814 }
 815 /**************************************************************************************************/
 816 int TrialSwap2::print_matrix(vector<vector<int> > &matrix, int nrows, int ncols)
 817 {
 818     try {
 819          m->mothurOut("matrix:");  m->mothurOutEndLine();
 820
 821         for (int i = 0; i < nrows; i++)
 822         {
 823             if (m->control_pressed) { return 0; }
 824             for (int j = 0; j < ncols; j++)
 825             {
 826                 m->mothurOut(toString(matrix[i][j]));
 827             }
 828             m->mothurOutEndLine();
 829         }
 830         return 0;
 831     }
 832     catch(exception& e) {
 833         m->errorOut(e, "TrialSwap2", "print_matrix");
 834         exit(1);
 835     }
 836 }
 837 /**************************************************************************************************/
 838
 839
 840
 841
 842
 843