//**********************************************************************************************************************
int CooccurrenceCommand::getCooccurrence(vector<SharedRAbundVector*>& thisLookUp, ofstream& out){
- try {
+ try {
int numOTUS = thisLookUp[0]->getNumBins();
vector< vector<int> > initmatrix; initmatrix.resize(thisLookUp.size());
vector< vector<int> > co_matrix; co_matrix.resize(thisLookUp[0]->getNumBins());
vector<int> rowtotal; rowtotal.resize(numOTUS, 0);
for (int i = 0; i < thisLookUp.size(); i++) {
- for (int j = 0; j < thisLookUp[i]->getNumBins(); j++) {
- if (m->control_pressed) { return 0; }
- int abund = thisLookUp[i]->getAbundance(j);
-
- if(abund > 0) {
- initmatrix[i][j] = 1;
+ for (int j = 0; j < thisLookUp[i]->getNumBins(); j++) {
+ if (m->control_pressed) { return 0; }
+ int abund = thisLookUp[i]->getAbundance(j);
+
+ if(abund > 0) {
+ initmatrix[i][j] = 1;
co_matrix[j][i] = 1;
rowtotal[j]++;
columntotal[i]++;
- }
- }
+ }
+ }
}
//nrows is ncols of inital matrix. All the functions need this value. They assume the transposition has already taken place and nrows and ncols refer to that matrix.
int nrows = numOTUS;//rows of inital matrix
int ncols = thisLookUp.size();//groups
double initscore = 0.0;
-
+
vector<int> columntotal; columntotal.resize(ncols, 0);
vector<int> rowtotal; rowtotal.resize(nrows, 0);
vector<double> stats;
double probabilityMatrix[ncols * nrows];
vector<vector<int> > nullmatrix(nrows, vector<int>(ncols, 0));
-
+
TrialSwap2 trial;
- n = accumulate( columntotal.begin(), columntotal.end(), 0 );
+ int n = accumulate( columntotal.begin(), columntotal.end(), 0 );
//============================================================
//generate a probability matrix. Only do this once.
float start = 0.0;
-
+
if (matrix == "sim1") {
for(int i=0;i<nrows;i++) {
for(int j=0;j<ncols;j++) {
- probabilityMatrix[ncols * i + j] = start + 1/double(nrows*ncols);
+ probabilityMatrix[ncols * i + j] = start + 1/double(nrows*ncols);
start = start + 1/double(nrows*ncols);
}
}
start = 0.0;
for(int j=0;j<ncols;j++) {
probabilityMatrix[ncols * i + j] = start + 1/double(ncols);
- start = start + 1/double(ncols);
+ start = start + 1/double(ncols);
}
}
}
-
+
else if (matrix == "sim3") {
for(int j=0;j<ncols;j++) {
start = 0.0;
for(int i=0;i<nrows;i++) {
- probabilityMatrix[ncols * i + j] = start + 1/double(nrows);
- start = start + 1/double(nrows);
+ probabilityMatrix[ncols * i + j] = start + 1/double(nrows);
+ start = start + 1/double(nrows);
}
}
}
-
+
else if (matrix == "sim4") {
for(int i=0;i<nrows;i++) {
start = 0.0;
for(int j=0;j<ncols;j++) {
- probabilityMatrix[ncols * i + j] = start + columntotal[j]/double(n);
- start = start + columntotal[j]/double(n);
+ probabilityMatrix[ncols * i + j] = start + columntotal[j]/double(n);
+ start = start + columntotal[j]/double(n);
}
}
}
-
+
else if (matrix == "sim5") {
for(int j=0;j<ncols;j++) {
start = 0.0;
for(int i=0;i<nrows;i++) {
- probabilityMatrix[ncols * i + j] = start + rowtotal[i]/double(n);
- start = start + rowtotal[i]/double(n);
+ probabilityMatrix[ncols * i + j] = start + rowtotal[i]/double(n);
+ start = start + rowtotal[i]/double(n);
}
}
}
-
+
else if (matrix == "sim6") {
for(int i=0;i<nrows;i++) {
for(int j=0;j<ncols;j++) {
- probabilityMatrix[ncols * i + j] = start + columntotal[j]/double(n*nrows);
- start = start + columntotal[j]/double(n*nrows);
+ probabilityMatrix[ncols * i + j] = start + columntotal[j]/double(n*nrows);
+ start = start + columntotal[j]/double(n*nrows);
}
}
}
-
-
+
+
else if (matrix == "sim7") {
for(int i=0;i<nrows;i++) {
for(int j=0;j<ncols;j++) {
- probabilityMatrix[ncols * i + j] = start + rowtotal[i]/double(n*ncols);
- start = start + rowtotal[i]/double(n*ncols);
+ probabilityMatrix[ncols * i + j] = start + rowtotal[i]/double(n*ncols);
+ start = start + rowtotal[i]/double(n*ncols);
}
}
}
-
+
else if (matrix == "sim8") {
for(int i=0;i<nrows;i++) {
for(int j=0;j<ncols;j++) {
- probabilityMatrix[ncols * i + j] = start + (rowtotal[i]*columntotal[j])/double(n*n);
- start = start + (rowtotal[i]*columntotal[j])/double(n*n);
+ probabilityMatrix[ncols * i + j] = start + (rowtotal[i]*columntotal[j])/double(n*n);
+ start = start + (rowtotal[i]*columntotal[j])/double(n*n);
}
}
}
-
+ else if (matrix == "sim9") { }
else {
- if(sim != 9) {
- m->mothurOut("[ERROR]: No model selected! \n");
- m->control_pressed = true;
- }
+ m->mothurOut("[ERROR]: No model selected! \n");
+ m->control_pressed = true;
}
-
-
- if (metric == "cscore") { initscore = trial.calc_c_score(initmatrix, rowtotal, ncols, nrows); }
- else if (metric == "checker") { initscore = trial.calc_checker(initmatrix, rowtotal, ncols, nrows); }
- else if (metric == "vratio") { initscore = trial.calc_vratio(nrows, ncols, rowtotal, columntotal); }
- else if (metric == "combo") { initscore = trial.calc_combo(nrows, ncols, initmatrix); }
- else { m->mothurOut("[ERROR]: No metric selected!\n"); m->control_pressed = true; return 1; }
+
+
+ if (metric == "cscore") { initscore = trial.calc_c_score(initmatrix, rowtotal, ncols, nrows); }
+ else if (metric == "checker") { initscore = trial.calc_checker(initmatrix, rowtotal, ncols, nrows); }
+ else if (metric == "vratio") { initscore = trial.calc_vratio(nrows, ncols, rowtotal, columntotal); }
+ else if (metric == "combo") { initscore = trial.calc_combo(nrows, ncols, initmatrix); }
+ else { m->mothurOut("[ERROR]: No metric selected!\n"); m->control_pressed = true; return 1; }
m->mothurOut("Initial c score: " + toString(initscore)); m->mothurOutEndLine();
double current;
double randnum;
int count;
-
+
//burn-in
- for(int i=0;i<10000;i++){
+ for(int i=0;i<10000;i++){
nullmatrix.clear();
//zero-fill the null matrix
nullmatrix.assign(nrows, vector<int>(ncols, 0));
-
+
if(matrix == "sim1" || matrix == "sim6" || matrix == "sim8" || matrix == "sim7") {
count = 0;
while(count < n) {
- nextnum:
+ nextnum:
previous = 0.0;
randnum = rand() / double(RAND_MAX);
for(int i=0;i<nrows;i++) {
- for(int j=0;j<ncols;j++) {
- current = probabilityMatrix[ncols * i + j];
+ for(int j=0;j<ncols;j++) {
+ current = probabilityMatrix[ncols * i + j];
if(randnum <= current && randnum > previous) {
nullmatrix[i][j] = 1;
count++;
if (count > n) break;
else
- goto nextnum;
+ goto nextnum;
}
previous = current;
}
}
}
}
-
+
else if(matrix == "sim2" || matrix == "sim4") {
- for(int i=0;i<nrows;i++) {
- previous = 0.0;
- count = 0;
- while(count < rowtotal[i]) {
- randnum = rand() / double(RAND_MAX);
- for(int j=0;j<ncols;j++) {
- current = probabilityMatrix[ncols * i + j];
- if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
- nullmatrix[i][j] = 1;
- count++;
- previous = 0.0;
- break;
- }
- previous = current;
+ for(int i=0;i<nrows;i++) {
+ previous = 0.0;
+ count = 0;
+ while(count < rowtotal[i]) {
+ randnum = rand() / double(RAND_MAX);
+ for(int j=0;j<ncols;j++) {
+ current = probabilityMatrix[ncols * i + j];
+ if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
+ nullmatrix[i][j] = 1;
+ count++;
+ previous = 0.0;
+ break;
}
+ previous = current;
}
}
+ }
}
-
+
else if(matrix == "sim3" || matrix == "sim5") {
//columns
for(int j=0;j<ncols;j++) {
}
- //populate null matrix from probability matrix, do this a lot.
- for(int i=0;i<runs;i++){
+ //populate null matrix from probability matrix, do this a lot.
+ for(int i=0;i<runs;i++){
nullmatrix.clear();
//zero-fill the null matrix
nullmatrix.assign(nrows, vector<int>(ncols, 0));
-
+
if(matrix == "sim1" || matrix == "sim6" || matrix == "sim8" || matrix == "sim7") {
count = 0;
while(count < n) {
- nextnum:
+ nextnum2:
previous = 0.0;
randnum = rand() / double(RAND_MAX);
for(int i=0;i<nrows;i++) {
- for(int j=0;j<ncols;j++) {
- current = probabilityMatrix[ncols * i + j];
+ for(int j=0;j<ncols;j++) {
+ current = probabilityMatrix[ncols * i + j];
if(randnum <= current && randnum > previous) {
nullmatrix[i][j] = 1;
count++;
if (count > n) break;
else
- goto nextnum;
+ goto nextnum2;
}
previous = current;
}
}
}
}
-
+
else if(matrix == "sim2" || matrix == "sim4") {
- for(int i=0;i<nrows;i++) {
- previous = 0.0;
- count = 0;
- while(count < rowtotal[i]) {
- randnum = rand() / double(RAND_MAX);
- for(int j=0;j<ncols;j++) {
- current = probabilityMatrix[ncols * i + j];
- if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
- nullmatrix[i][j] = 1;
- count++;
- previous = 0.0;
- break;
- }
- previous = current;
+ for(int i=0;i<nrows;i++) {
+ previous = 0.0;
+ count = 0;
+ while(count < rowtotal[i]) {
+ randnum = rand() / double(RAND_MAX);
+ for(int j=0;j<ncols;j++) {
+ current = probabilityMatrix[ncols * i + j];
+ if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
+ nullmatrix[i][j] = 1;
+ count++;
+ previous = 0.0;
+ break;
}
+ previous = current;
}
}
+ }
}
-
+
else if(matrix == "sim3" || matrix == "sim5") {
//columns
for(int j=0;j<ncols;j++) {
trial.swap_checkerboards (initmatrix, rowtotal, columntotal, ncols, nrows);
}
else {
- cout << "[ERROR]: No null model selected!\n" << endl;
+ m->mothurOut("[ERROR]: No null model selected!\n\n"); m->control_pressed = true;
return 1;
}
-
+
//run metric on null matrix and add score to the stats vector
if (metric == "cscore"){
- stats.push_back(trial.calc_c_score(nullmatrix, rowtotal, ncols, nrows));
+ stats.push_back(trial.calc_c_score(nullmatrix, rowtotal, ncols, nrows));
}
else if (metric == "checker") {
stats.push_back(trial.calc_checker(nullmatrix, rowtotal, ncols, nrows));
stats.push_back(trial.calc_combo(nrows, ncols, nullmatrix));
}
else {
- cout << "[ERROR]: No metric selected!\n" << endl;
+ m->mothurOut("[ERROR]: No metric selected!\n\n"); m->control_pressed = true;
return 1;
}
-
+
}
-
-
+
+
double total = 0.0;
- for (int i=0; i<stats.size();i++) { total+=stats[i]; }
+ for (int i=0; i<stats.size();i++) { total+=stats[i]; }
- double nullMean = double (total/(double)stats.size());
+ double nullMean = double (total/(double)stats.size());
m->mothurOutEndLine(); m->mothurOut("average metric score: " + toString(nullMean)); m->mothurOutEndLine();
double pvalue = 0.0;
- if (metric == "cscore" || metric == "checker") { pvalue = trial.calc_pvalue_greaterthan (stats, initscore); }
- else{ pvalue = trial.calc_pvalue_lessthan (stats, initscore); }
+ if (metric == "cscore" || metric == "checker") { pvalue = trial.calc_pvalue_greaterthan (stats, initscore); }
+ else{ pvalue = trial.calc_pvalue_lessthan (stats, initscore); }
m->mothurOut("pvalue: " + toString(pvalue)); m->mothurOutEndLine();
out << metric << '\t' << thisLookUp[0]->getLabel() << '\t' << nullMean << '\t' << pvalue << endl;
return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "CooccurrenceCommand", "Cooccurrence");
- exit(1);
- }
+ }
+ catch(exception& e) {
+ m->errorOut(e, "CooccurrenceCommand", "Cooccurrence");
+ exit(1);
+ }
}
//**********************************************************************************************************************
//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::sim1(vector<vector<int> > &co_matrix){
- try {
- vector<int> randRow;
- vector<vector<int> > tmpmatrix;
- int nrows = co_matrix.size();
- int ncols = co_matrix[0].size();
-
- //clear co_matrix
- // for(i=0;i<nrows;i++)
- // {
- // co_matrix.clear();
- // }
-
- //cout << "building matrix" << endl;
- for(int i=0;i<nrows;i++){
- if (m->control_pressed) { break; }
-
- for(int j=0;j<ncols;j++){
- double randNum = rand() / double(RAND_MAX);
- //cout << randNum << endl;
-
- if(randNum > 0.5) {
- randRow.push_back(1);
- }else{
- randRow.push_back(0);
- }
- }
- tmpmatrix.push_back(randRow);
- randRow.clear();
- //cout << endl;
- }
- co_matrix = tmpmatrix;
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "TrialSwap2", "sim1");
- exit(1);
- }
-}
-/**************************************************************************************************/
-/*
- *row sums fixed, columns equiprobable
- */
-void TrialSwap2::sim2(vector<vector<int> > &co_matrix)
-{
- try {
-
- for(int i=0;i<co_matrix.size();i++)
- {
- if (m->control_pressed) { break; }
- random_shuffle( co_matrix[i].begin(), co_matrix[i].end() );
- }
- }
- catch(exception& e) {
- m->errorOut(e, "TrialSwap2", "sim2");
- exit(1);
- }
-}
-/**************************************************************************************************/
-int TrialSwap2::sim2plus(vector<int> rowtotal, vector<vector<int> > &co_matrix)
-{
- try {
- int nrows = co_matrix.size();
- int ncols = co_matrix[0].size();
- double cellprob = 1.0/ncols;
- vector<double> cellprobvec;
- vector<int> tmprow;
- vector<vector<int> > tmpmatrix;
- //double randNum;
-
- double start = 0.0;
-
- for(int i=0; i<ncols; i++)
- {
- if (m->control_pressed) { return 0; }
- cellprobvec.push_back(start + cellprob);
- start = cellprobvec[i];
- }
-
- for(int i=0; i<nrows; i++)
- {
- tmprow.assign(ncols, 0);
-
- while( accumulate( tmprow.begin(), tmprow.end(), 0 ) < rowtotal[i])
- {
- if (m->control_pressed) { return 0; }
- double randNum = rand() / double(RAND_MAX);
- //cout << randNum << endl;
- if(randNum <= cellprobvec[0])
- {
- tmprow[0] = 1;
- continue;
- }
- for(int j=1;j<ncols;j++)
- {
- //cout << range[j] << endl;
- if(randNum <= cellprobvec[j] && randNum > cellprobvec[j-1] && tmprow[j] != 1)
- {
- tmprow[j] = 1;
- }
- }
- }
- tmpmatrix.push_back(tmprow);
- tmprow.clear();
- }
- co_matrix = tmpmatrix;
- tmpmatrix.clear();
- cellprobvec.clear();
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "TrialSwap2", "sim2plus");
- exit(1);
- }
-}
-/**************************************************************************************************/
-/*
- * same as sim2 but using initmatrix which is the initial co-occurrence matrix before transposition
- * may have to be changed depending on what matrix 'seed' is used. One way to use is to transpose
- * every null matrix before using an index and use the random matrix as a seed for the next null.
- */
-/**************************************************************************************************/
-void TrialSwap2::sim3(vector<vector<int> > &initmatrix)
-{
- try {
- for(int i=0;i<initmatrix.size();i++)
- {
- if (m->control_pressed) { break; }
- random_shuffle( initmatrix[i].begin(), initmatrix[i].end() );
- }
-
- }
- catch(exception& e) {
- m->errorOut(e, "TrialSwap2", "sim3");
- exit(1);
- }
-}
-/**************************************************************************************************/
-/*
- *
- *
- *
- */
-/**************************************************************************************************/
-int TrialSwap2::sim4(vector<int> columntotal, vector<int> rowtotal, vector<vector<int> > &co_matrix)
-{
- try {
- vector<double> colProb;
- vector<int> tmprow;//(ncols, 7);
- vector<vector<int> > tmpmatrix;
- vector<double> range;
- vector<double> randNums;
- int ncols = columntotal.size();
- int nrows = rowtotal.size();
- tmprow.clear();
-
- double colSum = accumulate( columntotal.begin(), columntotal.end(), 0 );
- //cout << "col sum: " << colSum << endl;
- for(int i=0;i<ncols;i++)
- {
- if (m->control_pressed) { return 0; }
- colProb.push_back(columntotal[i]/colSum);
- }
-
- double start = 0.0;
-
- for(int i=0;i<ncols;i++)
- {
- if (m->control_pressed) { return 0; }
- range.push_back(start + colProb[i]);
- start = range[i];
- }
-
- for(int i=0;i<nrows;i++)
- {
- tmprow.assign(ncols, 0);
- if (m->control_pressed) { return 0; }
-
- while ( accumulate( tmprow.begin(), tmprow.end(), 0 ) < rowtotal[i])
- {
- if (m->control_pressed) { return 0; }
-
- double randNum = rand() / double(RAND_MAX);
- if(randNum <= range[0])
- {
- tmprow[0] = 1;
- continue;
- }
- for(int j=1;j<ncols;j++)
- {
- if(randNum <= range[j] && randNum > range[j-1] && tmprow[j] != 1)
- {
- tmprow[j] = 1;
- }
-
- }
- }
- tmpmatrix.push_back(tmprow);
- tmprow.clear();
- }
-
- co_matrix = tmpmatrix;
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "TrialSwap2", "sim4");
- exit(1);
- }
-}
-/**************************************************************************************************/
-/*
- * inverse of sim4, MUST BE TRANSPOSED BEFORE CO-OCCURRENCE ANALYSIS
- *
- *
- */
-/**************************************************************************************************/
-int TrialSwap2::sim5(vector<int> initcolumntotal,vector<int> initrowtotal, vector<vector<int> > &initmatrix)
-{
- try {
- vector<double> colProb;
- vector<int> tmprow;//(ncols, 7);
- vector<vector<int> > tmpmatrix;
- vector<double> range;
- vector<double> randNums;
- int ncols = initcolumntotal.size();
- int nrows = initrowtotal.size();
-
- tmprow.clear();
-
- double colSum = accumulate( initcolumntotal.begin(), initcolumntotal.end(), 0 );
- //cout << "col sum: " << colSum << endl;
- for(int i=0;i<ncols;i++)
- {
- if (m->control_pressed) { return 0; }
- colProb.push_back(initcolumntotal[i]/colSum);
- }
-
- double start = 0.0;
-
- for(int i=0;i<ncols;i++)
- {
- if (m->control_pressed) { return 0; }
- range.push_back(start + colProb[i]);
- start = range[i];
- }
-
- for(int i=0;i<nrows;i++)
- {
- tmprow.assign(ncols, 0);
- if (m->control_pressed) { return 0; }
-
- while ( accumulate( tmprow.begin(), tmprow.end(), 0 ) < initrowtotal[i])
- {
- if (m->control_pressed) { return 0; }
-
- double randNum = rand() / double(RAND_MAX);
- if(randNum <= range[0])
- {
- tmprow[0] = 1;
- continue;
- }
- for(int j=1;j<ncols;j++)
- {
- if(randNum <= range[j] && randNum > range[j-1] && tmprow[j] != 1)
- {
- tmprow[j] = 1;
- }
-
- }
- }
- tmpmatrix.push_back(tmprow);
- tmprow.clear();
- }
-
- initmatrix = tmpmatrix;
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "TrialSwap2", "sim5");
- exit(1);
- }
-}
-/**************************************************************************************************/
-/*
- *
- *
- *
- */
-/**************************************************************************************************/
-int TrialSwap2::sim6(vector<int> columntotal, vector<vector<int> > &co_matrix)
-{
- try {
- vector<vector<int> > tmpmatrix;
- vector<double> colProb;
- vector<int> tmprow;
- vector<double> range;
- int ncols = columntotal.size();
- int nrows = co_matrix.size();
-
- int colSum = accumulate( columntotal.begin(), columntotal.end(), 0 );
-
- for(int i=0;i<ncols;i++)
- {
- if (m->control_pressed) { return 0; }
- colProb.push_back(columntotal[i]/double (colSum));
- }
-
- double start = 0.0;
-
- for(int i=0;i<ncols;i++)
- {
- if (m->control_pressed) { return 0; }
- range.push_back(start + colProb[i]);
- start = range[i];
- }
-
- for(int i=0;i<nrows;i++)
- {
- if (m->control_pressed) { return 0; }
- tmprow.assign(ncols, 0);
- int tmprowtotal;
- tmprowtotal = (rand() / double (RAND_MAX)) * 10;
- while ( tmprowtotal > ncols) {
- if (m->control_pressed) { return 0; }
- tmprowtotal = (rand() / double (RAND_MAX)) * 10;
- }
- //cout << tmprowtotal << endl;
- //cout << accumulate( tmprow.begin(), tmprow.end(), 0 ) << endl;
-
- while ( accumulate( tmprow.begin(), tmprow.end(), 0 ) < tmprowtotal)
- {
- if (m->control_pressed) { return 0; }
- double randNum = rand() / double(RAND_MAX);
- //cout << randNum << endl;
- if(randNum <= range[0])
- {
- tmprow[0] = 1;
- continue;
- }
- for(int j=1;j<ncols;j++)
- {
- //cout << range[j] << endl;
- if(randNum <= range[j] && randNum > range[j-1] && tmprow[j] != 1)
- {
- tmprow[j] = 1;
- }
-
- }
-
-
- }
-
- tmpmatrix.push_back(tmprow);
- tmprow.clear();
- }
-
- co_matrix = tmpmatrix;
- tmpmatrix.clear();
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "TrialSwap2", "sim6");
- exit(1);
- }
-}
-/**************************************************************************************************/
-/*
- * MUST BE TRANSPOSED BEFORE CO-OCCURRENCE ANALYSIS
- *
- *
- */
-/**************************************************************************************************/
-int TrialSwap2::sim7(vector<int> initrowtotal, vector<vector<int> > &co_matrix)
-{
- try {
- vector<vector<double> > probmatrix;
- vector<vector<int> > tmpmatrix;
- vector<double> colProb;
- vector<double> probrow;
- vector<int> tmprow;
- vector<double> range;
- double nc;
- int ncols = co_matrix[0].size(); int nrows = co_matrix.size();
-
- tmpmatrix.assign(nrows, vector<int>(ncols, 0.));
-
- int rowsum = accumulate( initrowtotal.begin(), initrowtotal.end(), 0 );
-
- nc = rowsum * ncols;
- //cout << nc << endl;
-
- //assign null matrix based on probabilities
-
- double start = 0.0; // don't reset start -- probs should be from 0-1 thoughout the entire matrix
-
- for(int i=0;i<nrows;i++)
- {
- if (m->control_pressed) { return 0; }
- //cout << initrowtotal[i]/double(nc) << endl;
- double cellprob = initrowtotal[i]/double(nc);
- //cout << cellprob << endl;
- for(int j=0;j<ncols;j++)
- {
-
- probrow.push_back(start + cellprob);
- //cout << probrow[j] << endl;
- //cout << start << endl;
- start = start + cellprob;
- }
- probmatrix.push_back(probrow);
- probrow.clear();
- }
-
-
- //while(tmprowsum < rowsum)
- //for(int k=0;k<rowsum;k++)
- int k = 0;
- while(k < rowsum)
- {
- if (m->control_pressed) { return 0; }
- done:
- //cout << k << endl;
- //tmprowsum = accumulate( tmprowtotal.begin(), tmprowtotal.end(), 0 );
- double randNum = rand() / double(RAND_MAX);
- //cout << randNum << "+" << endl;
- //special case for the first entry
- if(randNum <= probmatrix[0][0] && tmpmatrix[0][0] != 1)
- {
- tmpmatrix[0][0] = 1;
- k++;
- //cout << k << endl;
- continue;
- }
-
-
- for(int i=0;i<nrows;i++)
- {
- if (m->control_pressed) { return 0; }
- for(int j=0;j<ncols;j++)
- {
- //cout << probmatrix[i][j] << endl;
- if(randNum <= probmatrix[i][j] && randNum > probmatrix[i][j-1] && tmpmatrix[i][j] != 1)
- {
- tmpmatrix[i][j] = 1;
- k++;
- //cout << k << endl;
- goto done;
- }
- //else
- //k = k-1;
- }
-
- }
-
- }
-
- co_matrix = tmpmatrix;
- return 0;
- //build probibility matrix
- /* for(int i=0;i<nrows;i++)
- {
- for(int j=0;j<ncols;j++)
- {
- probrow.push_back(rowtotal[i]/nc);
- }
- probmatrix.pushback(probrow);
- probrow.clear;
- }
- */
-
- /* int colSum = accumulate( initcolumntotal.begin(), initcolumntotal.end(), 0 );
-
- for(int i=0;i<ncols;i++)
- {
- colProb.push_back(initcolumntotal[i]/double (colSum));
- }
-
- double start = 0.0;
-
- for(int i=0;i<ncols;i++)
- {
- range.push_back(start + colProb[i]);
- start = range[i];
- }
-
- for(int i=0;i<nrows;i++)
- {
- tmprow.assign(ncols, 0);
- int tmprowtotal;
- tmprowtotal = (rand() / double (RAND_MAX)) * 10;
- while ( tmprowtotal > ncols)
- tmprowtotal = (rand() / double (RAND_MAX)) * 10;
- //cout << tmprowtotal << endl;
- //cout << accumulate( tmprow.begin(), tmprow.end(), 0 ) << endl;
-
- while ( accumulate( tmprow.begin(), tmprow.end(), 0 ) < tmprowtotal)
- {
- double randNum = rand() / double(RAND_MAX);
- //cout << randNum << endl;
- if(randNum <= range[0])
- {
- tmprow[0] = 1;
- continue;
- }
- for(int j=1;j<ncols;j++)
- {
- //cout << range[j] << endl;
- if(randNum <= range[j] && randNum > range[j-1] && tmprow[j] != 1)
- {
- tmprow[j] = 1;
- }
- }
- }
-
- tmpmatrix.push_back(tmprow);
- tmprow.clear();
- }
-
-/**************************************************************************************************/
-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++;
+ }
}
}
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)
//need to transpose so we can compare rows (row-major order)
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)
{
try {
- int ncols = co_matrix[0].size(); int nrows = co_matrix.size();
+ 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; }}
-
+
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));
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 / commitdiff