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 double TrialSwap2::calc_c_score (vector<vector<int> > &co_matrix,vector<int>  rowtotal)
 648 {
 649     try {
 650         double cscore = 0.0;
 651         double maxD;
 652         double D;
 653         double normcscore = 0.0;
 654         int nonzeros = 0;
 655         int ncols = co_matrix[0].size(); int nrows = rowtotal.size();
 656         vector<vector<double> > s; s.resize(nrows);
 657         for (int i = 0; i < nrows; i++) { s[i].resize(nrows,0.0); }//only fill half the matrix
 658
 659
 660         for(int i=0;i<nrows-1;i++)
 661         {
 662
 663             for(int j=i+1;j<nrows;j++)
 664             {
 665                 if (m->control_pressed) { return 0; }
 666                 for(int k=0;k<ncols;k++)
 667                 {
 668                     if((co_matrix[i][k]==1)&&(co_matrix[j][k]==1)) //if both are 1s ie co-occurrence
 669                         s[i][j]++; //s counts co-occurrences
 670                 }
 671
 672                 //rowtotal[i] = A, rowtotal[j] = B, ncols = P, s[i][j] = J
 673                 cscore += (rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);///(nrows*(nrows-1)/2);
 674                 D = (rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);
 675
 676                 if(ncols < (rowtotal[i] + rowtotal[j]))
 677                 {
 678                     maxD = (ncols-rowtotal[i])*(ncols-rowtotal[j]);
 679                 }
 680                 else
 681                 {
 682                     maxD = rowtotal[i] * rowtotal[j];
 683                 }
 684
 685                 if(maxD != 0)
 686                 {
 687                     normcscore += D/maxD;
 688                     nonzeros++;
 689                 }
 690             }
 691         }
 692
 693         cscore = cscore/(double)(nrows*(nrows-1)/2);
 694         //cout << "normalized c score: " << normcscore/nonzeros << endl;
 695
 696         return cscore;
 697     }
 698         catch(exception& e) {
 699                 m->errorOut(e, "TrialSwap2", "calc_c_score");
 700                 exit(1);
 701         }
 702 }
 703 /**************************************************************************************************/
 704 int TrialSwap2::calc_checker (vector<vector<int> > &co_matrix, vector<int>  rowtotal)
 705 {
 706     try {
 707         int cunits=0;
 708         //int s[nrows][ncols];
 709         int ncols = co_matrix[0].size(); int nrows = rowtotal.size();
 710         vector<vector<int> > s; s.resize(nrows);
 711         for (int i = 0; i < nrows; i++) { s[i].resize(nrows,0); }//only fill half the matrix
 712
 713         for(int i=0;i<nrows-1;i++)
 714         {
 715             for(int j=i+1;j<nrows;j++)
 716             {
 717                 if (m->control_pressed) { return 0; }
 718                 //s[i][j]=0;
 719                 for(int k=0;k<ncols;k++)
 720                 {
 721                     //cout << s[i][j] << endl;
 722                     //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].
 723                     if((co_matrix[i][k]==1)&&(co_matrix[j][k]==1)) //if both are 1s ie co-occurrence
 724                         s[i][j]++; //s counts co-occurrences
 725
 726                 }
 727                 //cout << "rowtotal: " << rowtotal[i] << endl;
 728                 //cout << "co-occurrences: " << s[i][j] << endl;
 729                 //cunits+=(rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);
 730                 if (s[i][j] == 0)
 731                 {
 732                     cunits+=1;
 733                 }
 734                 //cunits+=s[i][j];
 735             }
 736         }
 737
 738         return cunits;
 739     }
 740         catch(exception& e) {
 741                 m->errorOut(e, "TrialSwap2", "calc_checker");
 742                 exit(1);
 743         }
 744 }
 745 /**************************************************************************************************/
 746 double TrialSwap2::calc_vratio (vector<int> rowtotal, vector<int> columntotal)
 747 {
 748     try {
 749         int nrows = rowtotal.size();
 750         int ncols = columntotal.size();
 751         int sumCol = accumulate(columntotal.begin(), columntotal.end(), 0 );
 752        // int sumRow = accumulate(rowtotal.begin(), rowtotal.end(), 0 );
 753
 754         double colAvg = (double) sumCol / (double) ncols;
 755  //       double rowAvg = (double) sumRow / (double) nrows;
 756
 757         double p = 0.0;
 758
 759  //       double totalRowVar = 0.0;
 760         double rowVar = 0.0;
 761         double colVar = 0.0;
 762
 763         for(int i=0;i<nrows;i++)
 764         {
 765             if (m->control_pressed) { return 0; }
 766             p = (double) rowtotal[i]/(double) ncols;
 767             rowVar += p * (1.0-p);
 768         }
 769
 770         for(int i=0;i<ncols;i++)
 771         {
 772             if (m->control_pressed) { return 0; }
 773             colVar += pow(((double) columntotal[i]-colAvg),2);
 774         }
 775
 776         colVar = (1.0/(double)ncols) * colVar;
 777
 778         return colVar/rowVar;
 779     }
 780     catch(exception& e) {
 781         m->errorOut(e, "TrialSwap2", "calc_vratio");
 782         exit(1);
 783     }
 784
 785 }
 786 /**************************************************************************************************/
 787 int TrialSwap2::calc_combo (vector<vector<int> > &initmatrix)
 788 {
 789     try {
 790         int initrows = initmatrix.size();
 791         int unique = 0;
 792         int match = 0;
 793         int matches = 0;
 794         for(int i=0;i<initrows;i++)
 795         {
 796             match = 0;
 797             for(int j=i+1;j<=initrows;j++)
 798             {
 799                 if (m->control_pressed) { return 0; }
 800                 if( (initmatrix[i] == initmatrix[j]))
 801                 {
 802                     match++;
 803                     matches++;
 804                     break;
 805                 }
 806             }
 807
 808             //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
 809             if (match == 0)
 810                 unique++;
 811         }
 812         return unique;
 813     }
 814     catch(exception& e) {
 815         m->errorOut(e, "TrialSwap2", "calc_combo");
 816         exit(1);
 817     }
 818 }
 819 /**************************************************************************************************/
 820 int TrialSwap2::swap_checkerboards (vector<vector<int> > &co_matrix)
 821 {
 822     try {
 823         int ncols = co_matrix[0].size(); int nrows = co_matrix.size();
 824         int i, j, k, l;
 825         i = m->getRandomIndex(nrows-1);
 826         while((j = m->getRandomIndex(nrows-1) ) == i ) {;if (m->control_pressed) { return 0; }}
 827         k = m->getRandomIndex(ncols-1);
 828         while((l = m->getRandomIndex(ncols-1)) == k ) {;if (m->control_pressed) { return 0; }}
 829
 830         //cout << co_matrix[i][k] << " " << co_matrix[j][l] << endl;
 831         //cout << co_matrix[i][l] << " " << co_matrix[j][k] << endl;
 832         //cout << co_matrix[i][l] << " " << co_matrix[j][k] << endl;
 833         //cout << co_matrix[i][l] << " " << co_matrix[j][k] << endl;
 834         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
 835         {
 836             co_matrix[i][k]=1-co_matrix[i][k];
 837             co_matrix[i][l]=1-co_matrix[i][l];
 838             co_matrix[j][k]=1-co_matrix[j][k];
 839             co_matrix[j][l]=1-co_matrix[j][l];
 840             //cout << "swapped!" << endl;
 841         }
 842         //cout << "i: " << i << " j: " << j << " k: " << " l: " << l << endl;
 843         return 0;
 844     }
 845     catch(exception& e) {
 846         m->errorOut(e, "TrialSwap2", "swap_checkerboards");
 847         exit(1);
 848     }
 849 }
 850 /**************************************************************************************************/
 851 double TrialSwap2::calc_pvalue_greaterthan (vector<double> scorevec, double initialscore)
 852 {
 853     try {
 854         int runs = scorevec.size();
 855         double p = 0.0;
 856         for( int i=0;i<runs;i++)
 857         {
 858             if (m->control_pressed) { return 0; }
 859             if(scorevec[i]>=initialscore)
 860                 p++;
 861         }
 862         return p/(double)runs;
 863     }
 864     catch(exception& e) {
 865         m->errorOut(e, "TrialSwap2", "calc_pvalue_greaterthan");
 866         exit(1);
 867     }
 868 }
 869 /**************************************************************************************************/
 870 double TrialSwap2::calc_pvalue_lessthan (vector<double> scorevec, double initialscore)
 871 {
 872     try {
 873         int runs = scorevec.size();
 874         double p = 0.0;
 875         for( int i=0;i<runs;i++)
 876         {
 877             if (m->control_pressed) { return 0; }
 878             if(scorevec[i]<=initialscore)
 879                 p++;
 880         }
 881         return p/(double)runs;
 882     }
 883     catch(exception& e) {
 884         m->errorOut(e, "TrialSwap2", "calc_pvalue_lessthan");
 885         exit(1);
 886     }
 887 }
 888 /**************************************************************************************************/
 889 double TrialSwap2::t_test (double initialscore, int runs, double nullMean, vector<double> scorevec)
 890 {
 891     try {
 892         double t;
 893         double sampleSD;
 894         double sum = 0;
 895
 896         for(int i=0;i<runs;i++)
 897         {
 898             if (m->control_pressed) { return 0; }
 899             sum += pow((scorevec[i] - nullMean),2);
 900             //cout << "scorevec[" << i << "]" << scorevec[i] << endl;
 901         }
 902
 903         m->mothurOut("nullMean: " + toString(nullMean)); m->mothurOutEndLine();
 904
 905         m->mothurOut("sum: " + toString(sum));  m->mothurOutEndLine();
 906
 907         sampleSD = sqrt( (1/runs) * sum );
 908
 909         m->mothurOut("samplSD: " + toString(sampleSD));  m->mothurOutEndLine();
 910
 911         t = (nullMean - initialscore) / (sampleSD / sqrt(runs));
 912
 913         return t;
 914     }
 915     catch(exception& e) {
 916         m->errorOut(e, "TrialSwap2", "t_test");
 917         exit(1);
 918     }
 919 }
 920 /**************************************************************************************************/
 921 int TrialSwap2::print_matrix(vector<vector<int> > &matrix, int nrows, int ncols)
 922 {
 923     try {
 924          m->mothurOut("matrix:");  m->mothurOutEndLine();
 925
 926         for (int i = 0; i < nrows; i++)
 927         {
 928             if (m->control_pressed) { return 0; }
 929             for (int j = 0; j < ncols; j++)
 930             {
 931                 m->mothurOut(toString(matrix[i][j]));
 932             }
 933             m->mothurOutEndLine();
 934         }
 935         return 0;
 936     }
 937     catch(exception& e) {
 938         m->errorOut(e, "TrialSwap2", "print_matrix");
 939         exit(1);
 940     }
 941 }
 942 /**************************************************************************************************/
 943 int TrialSwap2::transpose_matrix (vector<vector<int> > &initmatrix, vector<vector<int> > &co_matrix)//, int nrows, int nocols)
 944 {
 945     try {
 946         int ncols = initmatrix.size(); int nrows = initmatrix[0].size();
 947         int tmpnrows = nrows;
 948         //vector<vector<int> > tmpvec;
 949         vector<int> tmprow;
 950         if(!co_matrix.empty())
 951             co_matrix.clear();
 952         for (int i=0;i<nrows;i++)
 953         {
 954             if (m->control_pressed) { return 0; }
 955             for (int j=0;j<ncols;j++)
 956             {
 957                 tmprow.push_back(initmatrix[j][i]);
 958             }
 959             /*if (accumulate( tmprow.begin(), tmprow.end(), 0 ) == 0)
 960              {
 961              tmpnrows--;
 962              }
 963              else */
 964             co_matrix.push_back(tmprow);
 965             tmprow.clear();
 966         }
 967         nrows = tmpnrows;
 968         return 0;
 969     }
 970     catch(exception& e) {
 971         m->errorOut(e, "TrialSwap2", "transpose_matrix");
 972         exit(1);
 973     }
 974 }
 975 /**************************************************************************************************/
 976 int TrialSwap2::update_row_col_totals(vector<vector<int> > &co_matrix, vector<int> &rowtotal, vector<int> &columntotal)
 977 {
 978     try {
 979         //rowtotal.clear();
 980         //columntotal.clear();
 981         //generate (rowtotal.begin(), rowtotal.end(), 0);
 982         //generate (columntotal.begin(), columntotal.end(), 0);
 983         int nrows = co_matrix.size();
 984         int ncols = co_matrix[0].size();
 985         vector<int> tmpcolumntotal; tmpcolumntotal.resize(ncols, 0);
 986         vector<int> tmprowtotal; tmprowtotal.resize(nrows, 0);
 987
 988         int rowcount = 0;
 989
 990         for (int i = 0; i < nrows; i++)
 991         {
 992             if (m->control_pressed) { return 0; }
 993             for (int j = 0; j < ncols; j++)
 994             {
 995                 if (co_matrix[i][j] == 1)
 996                 {
 997                     rowcount++;
 998                     tmpcolumntotal[j]++;
 999                 }
1000             }
1001             tmprowtotal[i] = rowcount;
1002             rowcount = 0;
1003         }
1004         columntotal = tmpcolumntotal;
1005         rowtotal = tmprowtotal;
1006         /*cout << "rowtotal: ";
1007         for(int i = 0; i<nrows; i++) { cout << rowtotal[i]; }
1008         cout << "  ";
1009         cout << " coltotal: ";
1010         for(int i = 0; i<ncols; i++) { cout << columntotal[i]; }
1011         cout << endl;*/
1012         return 0;
1013     }
1014     catch(exception& e) {
1015         m->errorOut(e, "TrialSwap2", "update_row_col_totals");
1016         exit(1);
1017     }
1018 }
1019 /**************************************************************************************************/
1020
1021
1022
1023
1024