X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=trialSwap2.cpp;fp=trialSwap2.cpp;h=812077d16c98c5a3c865f197d1ba6f50fbce07a6;hb=41a0537df3cb09cb8c725eab126fcf9bf28ebd8d;hp=ea43de551d2a8df83abd5594754f17ebca2f6760;hpb=c669e63b8cb1193d683a1016562d35b455f76575;p=mothur.git

diff --git a/trialSwap2.cpp b/trialSwap2.cpp
index ea43de5..812077d 100644
--- a/trialSwap2.cpp
+++ b/trialSwap2.cpp
@@ -4,24 +4,528 @@
 //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){
+/**************************************************************************************************/
+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 {
-        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));
+        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", "intrand");
+		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)