//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.
-/**************************************************************************************************
-int TrialSwap2::intrand(int n){
- try {
- double z;
-
- z = (double)random() * (double)n / (double)RAND_MAX;
- if(z>=n)
- z=n-1;
- if(z<0)
- z=0;
- return((int)floor(z));
- }
- catch(exception& e) {
- m->errorOut(e, "TrialSwap2", "intrand");
- exit(1);
- }
-}
-/**************************************************************************************************/
-
-/**************************************************************************************************/
-double TrialSwap2::calc_c_score (vector<vector<int> > &co_matrix, vector<int> rowtotal, int ncols, int nrows)
+double TrialSwap2::calc_c_score (vector<vector<int> > &co_matrix, vector<int> rowtotal, int ncols, int nrows)
{
try {
double cscore = 0.0;
double D;
double normcscore = 0.0;
int nonzeros = 0;
- //int ncols = co_matrix[0].size(); int nrows = rowtotal.size();
+ //int ncols = co_matrix[0].size(); int nrows = rowtotal.size();
vector<vector<double> > s; s.resize(nrows);
for (int i = 0; i < nrows; i++) { s[i].resize(nrows,0.0); }//only fill half the matrix
-
+
for(int i=0;i<nrows-1;i++)
{
if(maxD != 0)
{
normcscore += D/maxD;
- nonzeros++;
- }
+ nonzeros++;
+ }
}
}
- cscore = cscore/(double)(nrows*(nrows-1)/2);
+ //cscore = cscore/(double)(nrows*(nrows-1)/2); //not normalized
//cout << "normalized c score: " << normcscore/nonzeros << endl;
-
+ cscore = normcscore/(double)nonzeros;
+
return cscore;
}
- catch(exception& e) {
- m->errorOut(e, "TrialSwap2", "calc_c_score");
- exit(1);
- }
+ catch(exception& e) {
+ m->errorOut(e, "TrialSwap2", "calc_c_score");
+ exit(1);
+ }
}
/**************************************************************************************************/
-int TrialSwap2::calc_checker (vector<vector<int> > &co_matrix, vector<int> rowtotal, int ncols, int nrows)
+int TrialSwap2::calc_checker (vector<vector<int> > &co_matrix, vector<int> rowtotal, int ncols, int nrows)
{
try {
int cunits=0;
//int s[nrows][ncols];
- //int ncols = co_matrix[0].size(); int nrows = rowtotal.size();
+ //int ncols = co_matrix[0].size(); int nrows = rowtotal.size();
vector<vector<int> > s; s.resize(nrows);
for (int i = 0; i < nrows; i++) { s[i].resize(nrows,0); }//only fill half the matrix
}
}
- return cunits;
+ return cunits;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "TrialSwap2", "calc_checker");
+ exit(1);
}
- catch(exception& e) {
- m->errorOut(e, "TrialSwap2", "calc_checker");
- exit(1);
- }
}
/**************************************************************************************************/
double TrialSwap2::calc_vratio (int nrows, int ncols, vector<int> rowtotal, vector<int> columntotal)
//int nrows = rowtotal.size();
//int ncols = columntotal.size();
int sumCol = accumulate(columntotal.begin(), columntotal.end(), 0 );
- // int sumRow = accumulate(rowtotal.begin(), rowtotal.end(), 0 );
+ // int sumRow = accumulate(rowtotal.begin(), rowtotal.end(), 0 );
double colAvg = (double) sumCol / (double) ncols;
- // double rowAvg = (double) sumRow / (double) nrows;
+ // double rowAvg = (double) sumRow / (double) nrows;
double p = 0.0;
- // double totalRowVar = 0.0;
+ // double totalRowVar = 0.0;
double rowVar = 0.0;
double colVar = 0.0;
if (m->control_pressed) { return 0; }
p = (double) rowtotal[i]/(double) ncols;
rowVar += p * (1.0-p);
- }
+ }
for(int i=0;i<ncols;i++)
{
m->errorOut(e, "TrialSwap2", "calc_vratio");
exit(1);
}
-
+
}
/**************************************************************************************************/
int TrialSwap2::calc_combo (int nrows, int ncols, vector<vector<int> > &nullmatrix)
{
try {
//need to transpose so we can compare rows (row-major order)
- int tmpnrows = nrows;
+ //int tmpnrows = nrows;
vector<vector<int> > tmpmatrix;
-
+
vector<int> tmprow;
if(!tmpmatrix.empty())
tmpmatrix.clear();
for (int i=0;i<ncols;i++)
- {
+ {
for (int j=0;j<nrows;j++)
{
tmprow.push_back(nullmatrix[j][i]);
}
-
+
tmpmatrix.push_back(tmprow);
tmprow.clear();
}
-
+
int unique = 0;
int match = 0;
for(int j=0;j<ncols;j++)
for(int i=j+1;i<=ncols;i++)
{
//comparing matrix rows
- if( (tmpmatrix[j] == tmpmatrix[i]))
+ if( (tmpmatrix[j] == tmpmatrix[i]))
{
match++;
break;
}
- }
-
+ }
+
//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
if (match == 0)
unique++;
- }
- return unique;
+ }
+ return unique;
}
catch(exception& e) {
m->errorOut(e, "TrialSwap2", "calc_combo");
exit(1);
}
-}
+}
/**************************************************************************************************/
-int TrialSwap2::swap_checkerboards (vector<vector<int> > &co_matrix)
+int TrialSwap2::swap_checkerboards (vector<vector<int> > &co_matrix, int ncols, int nrows)
{
try {
- int ncols = co_matrix[0].size(); int nrows = co_matrix.size();
- int i, j, k, l;
- i = m->getRandomIndex(nrows-1);
- while((j = m->getRandomIndex(nrows-1) ) == i ) {;if (m->control_pressed) { return 0; }}
- k = m->getRandomIndex(ncols-1);
- while((l = m->getRandomIndex(ncols-1)) == k ) {;if (m->control_pressed) { return 0; }}
+ //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.
+ for(int a=0;a<1000;a++){
+ int i, j, k, l;
+ i = m->getRandomIndex(nrows-1);
+ while((j = m->getRandomIndex(nrows-1) ) == i ) {;if (m->control_pressed) { return 0; }}
+ k = m->getRandomIndex(ncols-1);
+ while((l = m->getRandomIndex(ncols-1)) == k ) {;if (m->control_pressed) { return 0; }}
- 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
- {
- co_matrix[i][k]=1-co_matrix[i][k];
- co_matrix[i][l]=1-co_matrix[i][l];
- co_matrix[j][k]=1-co_matrix[j][k];
- co_matrix[j][l]=1-co_matrix[j][l];
+ 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
+ {
+ co_matrix[i][k]=1-co_matrix[i][k];
+ co_matrix[i][l]=1-co_matrix[i][l];
+ co_matrix[j][k]=1-co_matrix[j][k];
+ co_matrix[j][l]=1-co_matrix[j][l];
+ }
}
-
+
return 0;
}
catch(exception& e) {
m->mothurOut("nullMean: " + toString(nullMean)); m->mothurOutEndLine();
- m->mothurOut("sum: " + toString(sum)); m->mothurOutEndLine();
+ m->mothurOut("sum: " + toString(sum)); m->mothurOutEndLine();
sampleSD = sqrt( (1/runs) * sum );
- m->mothurOut("samplSD: " + toString(sampleSD)); m->mothurOutEndLine();
+ m->mothurOut("samplSD: " + toString(sampleSD)); m->mothurOutEndLine();
t = (nullMean - initialscore) / (sampleSD / sqrt(runs));
}
}
/**************************************************************************************************/
+double TrialSwap2::getSD (int runs, vector<double> scorevec, double nullMean)
+{
+ try{
+ double sum = 0;
+ for(int i=0;i<runs;i++)
+ {
+ if (m->control_pressed) { return 0; }
+ sum += pow((scorevec[i] - nullMean),2);
+ }
+ return sqrt( (1/double(runs)) * sum );
+ }
+ catch(exception& e) {
+ m->errorOut(e, "TrialSwap2", "getSD");
+ exit(1);
+ }
+}
+/**************************************************************************************************/
+double TrialSwap2::get_zscore (double sd, double nullMean, double initscore)
+{
+ try {
+ return (initscore - nullMean) / sd;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "TrialSwap2", "get_zscore");
+ exit(1);
+ }
+}
+/**************************************************************************************************/
int TrialSwap2::print_matrix(vector<vector<int> > &matrix, int nrows, int ncols)
{
try {
- m->mothurOut("matrix:"); m->mothurOutEndLine();
+ m->mothurOut("matrix:"); m->mothurOutEndLine();
for (int i = 0; i < nrows; i++)
{
if (m->control_pressed) { return 0; }
for (int j = 0; j < ncols; j++)
{
- m->mothurOut(toString(matrix[i][j]));
- }
+ m->mothurOut(toString(matrix[i][j]));
+ }
m->mothurOutEndLine();
}
return 0;
-
projects / mothur.git / blobdiff