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