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, int ncols, int nrows)
 196 {
 197     try {
 198         //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.
 199         for(int a=0;a<1000;a++){
 200             int i, j, k, l;
 201             i = m->getRandomIndex(nrows-1);
 202             while((j = m->getRandomIndex(nrows-1) ) == i ) {;if (m->control_pressed) { return 0; }}
 203             k = m->getRandomIndex(ncols-1);
 204             while((l = m->getRandomIndex(ncols-1)) == k ) {;if (m->control_pressed) { return 0; }}
 205
 206             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
 207             {
 208                 co_matrix[i][k]=1-co_matrix[i][k];
 209                 co_matrix[i][l]=1-co_matrix[i][l];
 210                 co_matrix[j][k]=1-co_matrix[j][k];
 211                 co_matrix[j][l]=1-co_matrix[j][l];
 212
 213             }
 214         }
 215
 216         return 0;
 217     }
 218     catch(exception& e) {
 219         m->errorOut(e, "TrialSwap2", "swap_checkerboards");
 220         exit(1);
 221     }
 222 }
 223 /**************************************************************************************************/
 224 double TrialSwap2::calc_pvalue_greaterthan (vector<double> scorevec, double initialscore)
 225 {
 226     try {
 227         int runs = scorevec.size();
 228         double p = 0.0;
 229         for( int i=0;i<runs;i++)
 230         {
 231             if (m->control_pressed) { return 0; }
 232             if(scorevec[i]>=initialscore)
 233                 p++;
 234         }
 235         return p/(double)runs;
 236     }
 237     catch(exception& e) {
 238         m->errorOut(e, "TrialSwap2", "calc_pvalue_greaterthan");
 239         exit(1);
 240     }
 241 }
 242 /**************************************************************************************************/
 243 double TrialSwap2::calc_pvalue_lessthan (vector<double> scorevec, double initialscore)
 244 {
 245     try {
 246         int runs = scorevec.size();
 247         double p = 0.0;
 248         for( int i=0;i<runs;i++)
 249         {
 250             if (m->control_pressed) { return 0; }
 251             if(scorevec[i]<=initialscore)
 252                 p++;
 253         }
 254         return p/(double)runs;
 255     }
 256     catch(exception& e) {
 257         m->errorOut(e, "TrialSwap2", "calc_pvalue_lessthan");
 258         exit(1);
 259     }
 260 }
 261 /**************************************************************************************************/
 262 double TrialSwap2::t_test (double initialscore, int runs, double nullMean, vector<double> scorevec)
 263 {
 264     try {
 265         double t;
 266         double sampleSD;
 267         double sum = 0;
 268
 269         for(int i=0;i<runs;i++)
 270         {
 271             if (m->control_pressed) { return 0; }
 272             sum += pow((scorevec[i] - nullMean),2);
 273             //cout << "scorevec[" << i << "]" << scorevec[i] << endl;
 274         }
 275
 276         m->mothurOut("nullMean: " + toString(nullMean)); m->mothurOutEndLine();
 277
 278         m->mothurOut("sum: " + toString(sum)); m->mothurOutEndLine();
 279
 280         sampleSD = sqrt( (1/runs) * sum );
 281
 282         m->mothurOut("samplSD: " + toString(sampleSD)); m->mothurOutEndLine();
 283
 284         t = (nullMean - initialscore) / (sampleSD / sqrt(runs));
 285
 286         return t;
 287     }
 288     catch(exception& e) {
 289         m->errorOut(e, "TrialSwap2", "t_test");
 290         exit(1);
 291     }
 292 }
 293 /**************************************************************************************************/
 294 double TrialSwap2::getSD (int runs, vector<double> scorevec, double nullMean)
 295 {
 296     try{
 297         double sum = 0;
 298         for(int i=0;i<runs;i++)
 299             {
 300                 if (m->control_pressed) { return 0; }
 301                 sum += pow((scorevec[i] - nullMean),2);
 302             }
 303         return sqrt( (1/double(runs)) * sum );
 304     }
 305     catch(exception& e) {
 306         m->errorOut(e, "TrialSwap2", "getSD");
 307         exit(1);
 308     }
 309 }
 310 /**************************************************************************************************/
 311 double TrialSwap2::get_zscore (double sd, double nullMean, double initscore)
 312 {
 313     try {
 314         return (initscore - nullMean) / sd;
 315     }
 316     catch(exception& e) {
 317         m->errorOut(e, "TrialSwap2", "get_zscore");
 318         exit(1);
 319     }
 320 }
 321 /**************************************************************************************************/
 322 int TrialSwap2::print_matrix(vector<vector<int> > &matrix, int nrows, int ncols)
 323 {
 324     try {
 325         m->mothurOut("matrix:"); m->mothurOutEndLine();
 326
 327         for (int i = 0; i < nrows; i++)
 328         {
 329             if (m->control_pressed) { return 0; }
 330             for (int j = 0; j < ncols; j++)
 331             {
 332                 m->mothurOut(toString(matrix[i][j]));
 333             }
 334             m->mothurOutEndLine();
 335         }
 336         return 0;
 337     }
 338     catch(exception& e) {
 339         m->errorOut(e, "TrialSwap2", "print_matrix");
 340         exit(1);
 341     }
 342 }
 343 /**************************************************************************************************/
 344
 345
 346
 347
 348