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
  26 /**************************************************************************************************/
  27 double TrialSwap2::calc_c_score (vector<vector<int> > &co_matrix, vector<int>  rowtotal, int ncols, int nrows)
  28 {
  29     try {
  30         double cscore = 0.0;
  31         double maxD;
  32         double D;
  33         double normcscore = 0.0;
  34         int nonzeros = 0;
  35         //int ncols = co_matrix[0].size(); int nrows = rowtotal.size();
  36         vector<vector<double> > s; s.resize(nrows);
  37         for (int i = 0; i < nrows; i++) { s[i].resize(nrows,0.0); }//only fill half the matrix
  38
  39
  40         for(int i=0;i<nrows-1;i++)
  41         {
  42
  43             for(int j=i+1;j<nrows;j++)
  44             {
  45                 if (m->control_pressed) { return 0; }
  46                 for(int k=0;k<ncols;k++)
  47                 {
  48                     if((co_matrix[i][k]==1)&&(co_matrix[j][k]==1)) //if both are 1s ie co-occurrence
  49                         s[i][j]++; //s counts co-occurrences
  50                 }
  51
  52                 //rowtotal[i] = A, rowtotal[j] = B, ncols = P, s[i][j] = J
  53                 cscore += (rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);///(nrows*(nrows-1)/2);
  54                 D = (rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);
  55
  56                 if(ncols < (rowtotal[i] + rowtotal[j]))
  57                 {
  58                     maxD = (ncols-rowtotal[i])*(ncols-rowtotal[j]);
  59                 }
  60                 else
  61                 {
  62                     maxD = rowtotal[i] * rowtotal[j];
  63                 }
  64
  65                 if(maxD != 0)
  66                 {
  67                     normcscore += D/maxD;
  68                     nonzeros++;
  69                 }
  70             }
  71         }
  72
  73         cscore = cscore/(double)(nrows*(nrows-1)/2);
  74         //cout << "normalized c score: " << normcscore/nonzeros << endl;
  75
  76         return cscore;
  77     }
  78         catch(exception& e) {
  79                 m->errorOut(e, "TrialSwap2", "calc_c_score");
  80                 exit(1);
  81         }
  82 }
  83 /**************************************************************************************************/
  84 int TrialSwap2::calc_checker (vector<vector<int> > &co_matrix, vector<int>  rowtotal, int ncols, int nrows)
  85 {
  86     try {
  87         int cunits=0;
  88         //int s[nrows][ncols];
  89         //int ncols = co_matrix[0].size(); int nrows = rowtotal.size();
  90         vector<vector<int> > s; s.resize(nrows);
  91         for (int i = 0; i < nrows; i++) { s[i].resize(nrows,0); }//only fill half the matrix
  92
  93         for(int i=0;i<nrows-1;i++)
  94         {
  95             for(int j=i+1;j<nrows;j++)
  96             {
  97                 if (m->control_pressed) { return 0; }
  98                 //s[i][j]=0;
  99                 for(int k=0;k<ncols;k++)
 100                 {
 101                     //cout << s[i][j] << endl;
 102                     //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].
 103                     if((co_matrix[i][k]==1)&&(co_matrix[j][k]==1)) //if both are 1s ie co-occurrence
 104                         s[i][j]++; //s counts co-occurrences
 105
 106                 }
 107                 //cout << "rowtotal: " << rowtotal[i] << endl;
 108                 //cout << "co-occurrences: " << s[i][j] << endl;
 109                 //cunits+=(rowtotal[i]-s[i][j])*(rowtotal[j]-s[i][j]);
 110                 if (s[i][j] == 0)
 111                 {
 112                     cunits+=1;
 113                 }
 114                 //cunits+=s[i][j];
 115             }
 116         }
 117
 118         return cunits;
 119     }
 120         catch(exception& e) {
 121                 m->errorOut(e, "TrialSwap2", "calc_checker");
 122                 exit(1);
 123         }
 124 }
 125 /**************************************************************************************************/
 126 double TrialSwap2::calc_vratio (int nrows, int ncols, vector<int> rowtotal, vector<int> columntotal)
 127 {
 128     try {
 129         //int nrows = rowtotal.size();
 130         //int ncols = columntotal.size();
 131         int sumCol = accumulate(columntotal.begin(), columntotal.end(), 0 );
 132        // int sumRow = accumulate(rowtotal.begin(), rowtotal.end(), 0 );
 133
 134         double colAvg = (double) sumCol / (double) ncols;
 135  //       double rowAvg = (double) sumRow / (double) nrows;
 136
 137         double p = 0.0;
 138
 139  //       double totalRowVar = 0.0;
 140         double rowVar = 0.0;
 141         double colVar = 0.0;
 142
 143         for(int i=0;i<nrows;i++)
 144         {
 145             if (m->control_pressed) { return 0; }
 146             p = (double) rowtotal[i]/(double) ncols;
 147             rowVar += p * (1.0-p);
 148         }
 149
 150         for(int i=0;i<ncols;i++)
 151         {
 152             if (m->control_pressed) { return 0; }
 153             colVar += pow(((double) columntotal[i]-colAvg),2);
 154         }
 155
 156         colVar = (1.0/(double)ncols) * colVar;
 157
 158         return colVar/rowVar;
 159     }
 160     catch(exception& e) {
 161         m->errorOut(e, "TrialSwap2", "calc_vratio");
 162         exit(1);
 163     }
 164
 165 }
 166 /**************************************************************************************************/
 167 int TrialSwap2::calc_combo (int nrows, int ncols, vector<vector<int> > &nullmatrix)
 168 {
 169     try {
 170         //need to transpose so we can compare rows (row-major order)
 171         int tmpnrows = nrows;
 172         vector<vector<int> > tmpmatrix;
 173
 174         vector<int> tmprow;
 175         if(!tmpmatrix.empty())
 176             tmpmatrix.clear();
 177         for (int i=0;i<ncols;i++)
 178         {
 179             for (int j=0;j<nrows;j++)
 180             {
 181                 tmprow.push_back(nullmatrix[j][i]);
 182             }
 183
 184             tmpmatrix.push_back(tmprow);
 185             tmprow.clear();
 186         }
 187
 188         int unique = 0;
 189         int match = 0;
 190         for(int j=0;j<ncols;j++)
 191         {
 192             match = 0;
 193             for(int i=j+1;i<=ncols;i++)
 194             {
 195                 //comparing matrix rows
 196                 if( (tmpmatrix[j] == tmpmatrix[i]))
 197                 {
 198                     match++;
 199                     break;
 200                 }
 201             }
 202
 203             //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
 204             if (match == 0)
 205                 unique++;
 206     }
 207     return unique;
 208     }
 209     catch(exception& e) {
 210         m->errorOut(e, "TrialSwap2", "calc_combo");
 211         exit(1);
 212     }
 213 }
 214 /**************************************************************************************************/
 215 int TrialSwap2::swap_checkerboards (vector<vector<int> > &co_matrix)
 216 {
 217     try {
 218         int ncols = co_matrix[0].size(); int nrows = co_matrix.size();
 219         int i, j, k, l;
 220         i = m->getRandomIndex(nrows-1);
 221         while((j = m->getRandomIndex(nrows-1) ) == i ) {;if (m->control_pressed) { return 0; }}
 222         k = m->getRandomIndex(ncols-1);
 223         while((l = m->getRandomIndex(ncols-1)) == k ) {;if (m->control_pressed) { return 0; }}
 224
 225         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
 226         {
 227             co_matrix[i][k]=1-co_matrix[i][k];
 228             co_matrix[i][l]=1-co_matrix[i][l];
 229             co_matrix[j][k]=1-co_matrix[j][k];
 230             co_matrix[j][l]=1-co_matrix[j][l];
 231
 232         }
 233
 234         return 0;
 235     }
 236     catch(exception& e) {
 237         m->errorOut(e, "TrialSwap2", "swap_checkerboards");
 238         exit(1);
 239     }
 240 }
 241 /**************************************************************************************************/
 242 double TrialSwap2::calc_pvalue_greaterthan (vector<double> scorevec, double initialscore)
 243 {
 244     try {
 245         int runs = scorevec.size();
 246         double p = 0.0;
 247         for( int i=0;i<runs;i++)
 248         {
 249             if (m->control_pressed) { return 0; }
 250             if(scorevec[i]>=initialscore)
 251                 p++;
 252         }
 253         return p/(double)runs;
 254     }
 255     catch(exception& e) {
 256         m->errorOut(e, "TrialSwap2", "calc_pvalue_greaterthan");
 257         exit(1);
 258     }
 259 }
 260 /**************************************************************************************************/
 261 double TrialSwap2::calc_pvalue_lessthan (vector<double> scorevec, double initialscore)
 262 {
 263     try {
 264         int runs = scorevec.size();
 265         double p = 0.0;
 266         for( int i=0;i<runs;i++)
 267         {
 268             if (m->control_pressed) { return 0; }
 269             if(scorevec[i]<=initialscore)
 270                 p++;
 271         }
 272         return p/(double)runs;
 273     }
 274     catch(exception& e) {
 275         m->errorOut(e, "TrialSwap2", "calc_pvalue_lessthan");
 276         exit(1);
 277     }
 278 }
 279 /**************************************************************************************************/
 280 double TrialSwap2::t_test (double initialscore, int runs, double nullMean, vector<double> scorevec)
 281 {
 282     try {
 283         double t;
 284         double sampleSD;
 285         double sum = 0;
 286
 287         for(int i=0;i<runs;i++)
 288         {
 289             if (m->control_pressed) { return 0; }
 290             sum += pow((scorevec[i] - nullMean),2);
 291             //cout << "scorevec[" << i << "]" << scorevec[i] << endl;
 292         }
 293
 294         m->mothurOut("nullMean: " + toString(nullMean)); m->mothurOutEndLine();
 295
 296         m->mothurOut("sum: " + toString(sum));  m->mothurOutEndLine();
 297
 298         sampleSD = sqrt( (1/runs) * sum );
 299
 300         m->mothurOut("samplSD: " + toString(sampleSD));  m->mothurOutEndLine();
 301
 302         t = (nullMean - initialscore) / (sampleSD / sqrt(runs));
 303
 304         return t;
 305     }
 306     catch(exception& e) {
 307         m->errorOut(e, "TrialSwap2", "t_test");
 308         exit(1);
 309     }
 310 }
 311 /**************************************************************************************************/
 312 int TrialSwap2::print_matrix(vector<vector<int> > &matrix, int nrows, int ncols)
 313 {
 314     try {
 315          m->mothurOut("matrix:");  m->mothurOutEndLine();
 316
 317         for (int i = 0; i < nrows; i++)
 318         {
 319             if (m->control_pressed) { return 0; }
 320             for (int j = 0; j < ncols; j++)
 321             {
 322                 m->mothurOut(toString(matrix[i][j]));
 323             }
 324             m->mothurOutEndLine();
 325         }
 326         return 0;
 327     }
 328     catch(exception& e) {
 329         m->errorOut(e, "TrialSwap2", "print_matrix");
 330         exit(1);
 331     }
 332 }
 333 /**************************************************************************************************/
 334
 335
 336
 337
 338
 339