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         initmatrix = tmpmatrix;
 531      */
 532     }
 533         catch(exception& e) {
 534                 m->errorOut(e, "TrialSwap2", "sim7");
 535                 exit(1);
 536         }
 537 }
 538 /**************************************************************************************************/
 539 /*
 540  *
 541  *
 542  *
 543  */
 544 /**************************************************************************************************/
 545 int TrialSwap2::sim8(vector<int> columntotal, vector<int> rowtotal, vector<vector<int> > &co_matrix)
 546 {
 547     try {
 548         double prob;
 549         double start = 0.0;
 550         int ncols = columntotal.size(); int nrows = rowtotal.size();
 551         double probarray[nrows * ncols];
 552         double randnum;
 553         int grandtotal;
 554         int total = 0;
 555
 556         //double colSum = accumulate( columntotal.begin(), columntotal.end(), 0 );
 557         double rowSum = accumulate( rowtotal.begin(), rowtotal.end(), 0 );
 558
 559         if (m->control_pressed) { return 0; }
 560
 561         //cout << "rowsum: " << rowSum << endl;
 562
 563         grandtotal = rowSum;
 564
 565         //create probability matrix with each site being between 0 and 1
 566         for (int i=0;i<nrows;i++) {
 567             if (m->control_pressed) { return 0; }
 568             for (int j=0;j<ncols;j++) {
 569                 prob = (rowtotal[i] * columntotal[j])/(rowSum*rowSum);
 570                 if (prob == 0.0)
 571                     probarray[ncols * i + j] = -1;
 572                 else
 573                     probarray[ncols * i + j] = start + prob;
 574                 //probmatrixrow.pushback(start + prob);
 575                 start += prob;
 576             }
 577         }
 578         //cout << "prbarray" << endl;
 579         //for(int i=0;i<(nrows*ncols);i++)
 580         //cout << probarray[i] << " ";
 581         //cout << endl;
 582
 583         //generate random muber between 0 and 1 and interate through probarray until found
 584         while (total < grandtotal)  {
 585             if (m->control_pressed) { return 0; }
 586             randnum = rand() / double(RAND_MAX);
 587             //cout << "rand num: " << randnum << endl;
 588             if((randnum <= probarray[0]) && (probarray[0] != 2) ) {
 589                 probarray[0] = 2;
 590                 total++;
 591                 continue;
 592             }
 593             for(int i=1;i<(nrows*ncols);i++) {
 594                 if (m->control_pressed) { return 0; }
 595                 if((randnum <= probarray[i]) && (randnum > probarray[i-1]) && (probarray[i] != 2) ) {
 596                     probarray[i] = 2;
 597                     total++;
 598                     break;
 599                 }
 600                 else
 601                     continue;
 602             }
 603         }
 604         //cout << "prbarray" << endl;
 605         //for(int i=0;i<(nrows*ncols);i++)
 606         //cout << probarray[i] << " ";
 607         //cout << endl;
 608         for(int i=0;i<nrows;i++) {
 609             if (m->control_pressed) { return 0; }
 610             for(int j=0;j<ncols;j++) {
 611                 if(probarray[ncols * i + j] == 2)
 612                     co_matrix[i][j] = 1;
 613                 else
 614                     co_matrix[i][j] = 0;
 615             }
 616         }
 617         return 0;
 618     }
 619         catch(exception& e) {
 620                 m->errorOut(e, "TrialSwap2", "sim8");
 621                 exit(1);
 622         }
 623 }
 624 /**************************************************************************************************/
 625 double TrialSwap2::calc_c_score (vector<vector<int> > &co_matrix,vector<int>  rowtotal)
 626 {
 627     try {
 628         double cscore = 0.0;
 629         double maxD;
 630         double D;
 631         double normcscore = 0.0;
 632         int nonzeros = 0;
 633         int ncols = co_matrix[0].size(); int nrows = rowtotal.size();
 634         vector<vector<double> > s; s.resize(nrows);
 635         for (int i = 0; i < nrows; i++) { s[i].resize(nrows,0.0); }//only fill half the matrix
 636
 637
 638         for(int i=0;i<nrows-1;i++)
 639         {
 640
 641             for(int j=i+1;j<nrows;j++)
 642             {
 643                 if (m->control_pressed) { return 0; }
 644                 for(int k=0;k<ncols;k++)
 645                 {
 646                     if((co_matrix[i][k]==1)&&(co_matrix[j][k]==1)) //if both are 1s ie co-occurrence
 647                         s[i][j]++; //s counts co-occurrences
 648                 }
 649
 650                 //rowtotal[i] = A, rowtotal[j] = B, ncols = P, s[i][j] = J
 651                 cscore += (rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);///(nrows*(nrows-1)/2);
 652                 D = (rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);
 653
 654                 if(ncols < (rowtotal[i] + rowtotal[j]))
 655                 {
 656                     maxD = (ncols-rowtotal[i])*(ncols-rowtotal[j]);
 657                 }
 658                 else
 659                 {
 660                     maxD = rowtotal[i] * rowtotal[j];
 661                 }
 662
 663                 if(maxD != 0)
 664                 {
 665                     normcscore += D/maxD;
 666                     nonzeros++;
 667                 }
 668             }
 669         }
 670
 671         cscore = cscore/(double)(nrows*(nrows-1)/2);
 672         //cout << "normalized c score: " << normcscore/nonzeros << endl;
 673
 674         return cscore;
 675     }
 676         catch(exception& e) {
 677                 m->errorOut(e, "TrialSwap2", "calc_c_score");
 678                 exit(1);
 679         }
 680 }
 681 /**************************************************************************************************/
 682 int TrialSwap2::calc_checker (vector<vector<int> > &co_matrix, vector<int>  rowtotal)
 683 {
 684     try {
 685         int cunits=0;
 686         //int s[nrows][ncols];
 687         int ncols = co_matrix[0].size(); int nrows = rowtotal.size();
 688         vector<vector<int> > s; s.resize(nrows);
 689         for (int i = 0; i < nrows; i++) { s[i].resize(nrows,0); }//only fill half the matrix
 690
 691         for(int i=0;i<nrows-1;i++)
 692         {
 693             for(int j=i+1;j<nrows;j++)
 694             {
 695                 if (m->control_pressed) { return 0; }
 696                 //s[i][j]=0;
 697                 for(int k=0;k<ncols;k++)
 698                 {
 699                     //cout << s[i][j] << endl;
 700                     //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].
 701                     if((co_matrix[i][k]==1)&&(co_matrix[j][k]==1)) //if both are 1s ie co-occurrence
 702                         s[i][j]++; //s counts co-occurrences
 703
 704                 }
 705                 //cout << "rowtotal: " << rowtotal[i] << endl;
 706                 //cout << "co-occurrences: " << s[i][j] << endl;
 707                 //cunits+=(rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);
 708                 if (s[i][j] == 0)
 709                 {
 710                     cunits+=1;
 711                 }
 712                 //cunits+=s[i][j];
 713             }
 714         }
 715
 716         return cunits;
 717     }
 718         catch(exception& e) {
 719                 m->errorOut(e, "TrialSwap2", "calc_checker");
 720                 exit(1);
 721         }
 722 }
 723 /**************************************************************************************************/
 724 double TrialSwap2::calc_vratio (vector<int> rowtotal, vector<int> columntotal)
 725 {
 726     try {
 727         int nrows = rowtotal.size();
 728         int ncols = columntotal.size();
 729         int sumCol = accumulate(columntotal.begin(), columntotal.end(), 0 );
 730        // int sumRow = accumulate(rowtotal.begin(), rowtotal.end(), 0 );
 731
 732         double colAvg = (double) sumCol / (double) ncols;
 733  //       double rowAvg = (double) sumRow / (double) nrows;
 734
 735         double p = 0.0;
 736
 737  //       double totalRowVar = 0.0;
 738         double rowVar = 0.0;
 739         double colVar = 0.0;
 740
 741         for(int i=0;i<nrows;i++)
 742         {
 743             if (m->control_pressed) { return 0; }
 744             p = (double) rowtotal[i]/(double) ncols;
 745             rowVar += p * (1.0-p);
 746         }
 747
 748         for(int i=0;i<ncols;i++)
 749         {
 750             if (m->control_pressed) { return 0; }
 751             colVar += pow(((double) columntotal[i]-colAvg),2);
 752         }
 753
 754         colVar = (1.0/(double)ncols) * colVar;
 755
 756         return colVar/rowVar;
 757     }
 758     catch(exception& e) {
 759         m->errorOut(e, "TrialSwap2", "calc_vratio");
 760         exit(1);
 761     }
 762
 763 }
 764 /**************************************************************************************************/
 765 int TrialSwap2::calc_combo (vector<vector<int> > &initmatrix)
 766 {
 767     try {
 768         int initrows = initmatrix.size();
 769         int unique = 0;
 770         int match = 0;
 771         int matches = 0;
 772         for(int i=0;i<initrows;i++)
 773         {
 774             match = 0;
 775             for(int j=i+1;j<=initrows;j++)
 776             {
 777                 if (m->control_pressed) { return 0; }
 778                 if( (initmatrix[i] == initmatrix[j]))
 779                 {
 780                     match++;
 781                     matches++;
 782                     break;
 783                 }
 784             }
 785
 786             //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
 787             if (match == 0)
 788                 unique++;
 789         }
 790         return unique;
 791     }
 792     catch(exception& e) {
 793         m->errorOut(e, "TrialSwap2", "calc_combo");
 794         exit(1);
 795     }
 796 }
 797 /**************************************************************************************************/
 798 int TrialSwap2::swap_checkerboards (vector<vector<int> > &co_matrix)
 799 {
 800     try {
 801         int ncols = co_matrix[0].size(); int nrows = co_matrix.size();
 802         int i, j, k, l;
 803         i = m->getRandomIndex(nrows-1);
 804         while((j = m->getRandomIndex(nrows-1) ) == i ) {;if (m->control_pressed) { return 0; }}
 805         k = m->getRandomIndex(ncols-1);
 806         while((l = m->getRandomIndex(ncols-1)) == k ) {;if (m->control_pressed) { return 0; }}
 807
 808         //cout << co_matrix[i][k] << " " << co_matrix[j][l] << endl;
 809         //cout << co_matrix[i][l] << " " << co_matrix[j][k] << endl;
 810         //cout << co_matrix[i][l] << " " << co_matrix[j][k] << endl;
 811         //cout << co_matrix[i][l] << " " << co_matrix[j][k] << endl;
 812         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
 813         {
 814             co_matrix[i][k]=1-co_matrix[i][k];
 815             co_matrix[i][l]=1-co_matrix[i][l];
 816             co_matrix[j][k]=1-co_matrix[j][k];
 817             co_matrix[j][l]=1-co_matrix[j][l];
 818             //cout << "swapped!" << endl;
 819         }
 820         //cout << "i: " << i << " j: " << j << " k: " << " l: " << l << endl;
 821         return 0;
 822     }
 823     catch(exception& e) {
 824         m->errorOut(e, "TrialSwap2", "swap_checkerboards");
 825         exit(1);
 826     }
 827 }
 828 /**************************************************************************************************/
 829 double TrialSwap2::calc_pvalue_greaterthan (vector<double> scorevec, double initialscore)
 830 {
 831     try {
 832         int runs = scorevec.size();
 833         double p = 0.0;
 834         for( int i=0;i<runs;i++)
 835         {
 836             if (m->control_pressed) { return 0; }
 837             if(scorevec[i]>=initialscore)
 838                 p++;
 839         }
 840         return p/(double)runs;
 841     }
 842     catch(exception& e) {
 843         m->errorOut(e, "TrialSwap2", "calc_pvalue_greaterthan");
 844         exit(1);
 845     }
 846 }
 847 /**************************************************************************************************/
 848 double TrialSwap2::calc_pvalue_lessthan (vector<double> scorevec, double initialscore)
 849 {
 850     try {
 851         int runs = scorevec.size();
 852         double p = 0.0;
 853         for( int i=0;i<runs;i++)
 854         {
 855             if (m->control_pressed) { return 0; }
 856             if(scorevec[i]<=initialscore)
 857                 p++;
 858         }
 859         return p/(double)runs;
 860     }
 861     catch(exception& e) {
 862         m->errorOut(e, "TrialSwap2", "calc_pvalue_lessthan");
 863         exit(1);
 864     }
 865 }
 866 /**************************************************************************************************/
 867 double TrialSwap2::t_test (double initialscore, int runs, double nullMean, vector<double> scorevec)
 868 {
 869     try {
 870         double t;
 871         double sampleSD;
 872         double sum = 0;
 873
 874         for(int i=0;i<runs;i++)
 875         {
 876             if (m->control_pressed) { return 0; }
 877             sum += pow((scorevec[i] - nullMean),2);
 878             //cout << "scorevec[" << i << "]" << scorevec[i] << endl;
 879         }
 880
 881         m->mothurOut("nullMean: " + toString(nullMean)); m->mothurOutEndLine();
 882
 883         m->mothurOut("sum: " + toString(sum));  m->mothurOutEndLine();
 884
 885         sampleSD = sqrt( (1/runs) * sum );
 886
 887         m->mothurOut("samplSD: " + toString(sampleSD));  m->mothurOutEndLine();
 888
 889         t = (nullMean - initialscore) / (sampleSD / sqrt(runs));
 890
 891         return t;
 892     }
 893     catch(exception& e) {
 894         m->errorOut(e, "TrialSwap2", "t_test");
 895         exit(1);
 896     }
 897 }
 898 /**************************************************************************************************/
 899 int TrialSwap2::print_matrix(vector<vector<int> > &matrix, int nrows, int ncols)
 900 {
 901     try {
 902          m->mothurOut("matrix:");  m->mothurOutEndLine();
 903
 904         for (int i = 0; i < nrows; i++)
 905         {
 906             if (m->control_pressed) { return 0; }
 907             for (int j = 0; j < ncols; j++)
 908             {
 909                 m->mothurOut(toString(matrix[i][j]));
 910             }
 911             m->mothurOutEndLine();
 912         }
 913         return 0;
 914     }
 915     catch(exception& e) {
 916         m->errorOut(e, "TrialSwap2", "print_matrix");
 917         exit(1);
 918     }
 919 }
 920 /**************************************************************************************************/
 921 int TrialSwap2::transpose_matrix (vector<vector<int> > &initmatrix, vector<vector<int> > &co_matrix)//, int nrows, int nocols)
 922 {
 923     try {
 924         int ncols = initmatrix.size(); int nrows = initmatrix[0].size();
 925         int tmpnrows = nrows;
 926         //vector<vector<int> > tmpvec;
 927         vector<int> tmprow;
 928         if(!co_matrix.empty())
 929             co_matrix.clear();
 930         for (int i=0;i<nrows;i++)
 931         {
 932             if (m->control_pressed) { return 0; }
 933             for (int j=0;j<ncols;j++)
 934             {
 935                 tmprow.push_back(initmatrix[j][i]);
 936             }
 937             /*if (accumulate( tmprow.begin(), tmprow.end(), 0 ) == 0)
 938              {
 939              tmpnrows--;
 940              }
 941              else */
 942             co_matrix.push_back(tmprow);
 943             tmprow.clear();
 944         }
 945         nrows = tmpnrows;
 946         return 0;
 947     }
 948     catch(exception& e) {
 949         m->errorOut(e, "TrialSwap2", "transpose_matrix");
 950         exit(1);
 951     }
 952 }
 953 /**************************************************************************************************/
 954 int TrialSwap2::update_row_col_totals(vector<vector<int> > &co_matrix, vector<int> &rowtotal, vector<int> &columntotal)
 955 {
 956     try {
 957         //rowtotal.clear();
 958         //columntotal.clear();
 959         //generate (rowtotal.begin(), rowtotal.end(), 0);
 960         //generate (columntotal.begin(), columntotal.end(), 0);
 961         int nrows = co_matrix.size();
 962         int ncols = co_matrix[0].size();
 963         vector<int> tmpcolumntotal; tmpcolumntotal.resize(ncols, 0);
 964         vector<int> tmprowtotal; tmprowtotal.resize(nrows, 0);
 965
 966         int rowcount = 0;
 967
 968         for (int i = 0; i < nrows; i++)
 969         {
 970             if (m->control_pressed) { return 0; }
 971             for (int j = 0; j < ncols; j++)
 972             {
 973                 if (co_matrix[i][j] == 1)
 974                 {
 975                     rowcount++;
 976                     tmpcolumntotal[j]++;
 977                 }
 978             }
 979             tmprowtotal[i] = rowcount;
 980             rowcount = 0;
 981         }
 982         columntotal = tmpcolumntotal;
 983         rowtotal = tmprowtotal;
 984         /*cout << "rowtotal: ";
 985         for(int i = 0; i<nrows; i++) { cout << rowtotal[i]; }
 986         cout << "  ";
 987         cout << " coltotal: ";
 988         for(int i = 0; i<ncols; i++) { cout << columntotal[i]; }
 989         cout << endl;*/
 990         return 0;
 991     }
 992     catch(exception& e) {
 993         m->errorOut(e, "TrialSwap2", "update_row_col_totals");
 994         exit(1);
 995     }
 996 }
 997 /**************************************************************************************************/
 998
 999
1000
1001
1002