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