+ #endif
+ }
+
+ if (subsample && (thisIter != 0)) {
+ if((thisIter) % 100 == 0){ m->mothurOut(toString(thisIter)); m->mothurOutEndLine(); }
+ calcDistsTotals.push_back(calcDists);
+ for (int i = 0; i < calcDists.size(); i++) {
+ for (int j = 0; j < calcDists[i].size(); j++) {
+ if (m->debug) { m->mothurOut("[DEBUG]: Results: iter = " + toString(thisIter) + ", " + thisLookup[calcDists[i][j].seq1]->getGroup() + " - " + thisLookup[calcDists[i][j].seq2]->getGroup() + " distance = " + toString(calcDists[i][j].dist) + ".\n"); }
+ }
+ }
+ //clean up memory
+ for (int i = 0; i < thisItersLookup.size(); i++) { delete thisItersLookup[i]; }
+ thisItersLookup.clear();
+ }else { //print results for whole dataset
+ for (int i = 0; i < calcDists.size(); i++) {
+ if (m->control_pressed) { break; }
+
+ //initialize matrix
+ vector< vector<double> > matrix; //square matrix to represent the distance
+ matrix.resize(thisLookup.size());
+ for (int k = 0; k < thisLookup.size(); k++) { matrix[k].resize(thisLookup.size(), 0.0); }
+
+ for (int j = 0; j < calcDists[i].size(); j++) {
+ int row = calcDists[i][j].seq1;
+ int column = calcDists[i][j].seq2;
+ double dist = calcDists[i][j].dist;
+
+ matrix[row][column] = dist;
+ matrix[column][row] = dist;
+ }
+
+ variables["[outputtag]"] = output;
+ variables["[calc]"] = matrixCalculators[i]->getName();
+ string distFileName = getOutputFileName("phylip",variables);
+ outputNames.push_back(distFileName); outputTypes["phylip"].push_back(distFileName);
+
+ ofstream outDist;
+ m->openOutputFile(distFileName, outDist);
+ outDist.setf(ios::fixed, ios::floatfield); outDist.setf(ios::showpoint);
+
+ printSims(outDist, matrix);
+
+ outDist.close();
+ }
+ }
+ for (int i = 0; i < calcDists.size(); i++) { calcDists[i].clear(); }
+ }
+
+ if (iters != 0) {
+ //we need to find the average distance and standard deviation for each groups distance
+
+ vector< vector<seqDist> > calcAverages; calcAverages.resize(matrixCalculators.size());
+ for (int i = 0; i < calcAverages.size(); i++) { //initialize sums to zero.
+ calcAverages[i].resize(calcDistsTotals[0][i].size());
+
+ for (int j = 0; j < calcAverages[i].size(); j++) {
+ calcAverages[i][j].seq1 = calcDistsTotals[0][i][j].seq1;
+ calcAverages[i][j].seq2 = calcDistsTotals[0][i][j].seq2;
+ calcAverages[i][j].dist = 0.0;
+ }
+ }
+ if (mode == "average") {
+ for (int thisIter = 0; thisIter < iters; thisIter++) { //sum all groups dists for each calculator
+ for (int i = 0; i < calcAverages.size(); i++) { //initialize sums to zero.
+ for (int j = 0; j < calcAverages[i].size(); j++) {
+ calcAverages[i][j].dist += calcDistsTotals[thisIter][i][j].dist;
+ if (m->debug) { m->mothurOut("[DEBUG]: Totaling for average calc: iter = " + toString(thisIter) + ", " + thisLookup[calcDistsTotals[thisIter][i][j].seq1]->getGroup() + " - " + thisLookup[calcDistsTotals[thisIter][i][j].seq2]->getGroup() + " distance = " + toString(calcDistsTotals[thisIter][i][j].dist) + ". New total = " + toString(calcAverages[i][j].dist) + ".\n"); }
+ }
+ }
+ }
+
+ for (int i = 0; i < calcAverages.size(); i++) { //finds average.
+ for (int j = 0; j < calcAverages[i].size(); j++) {
+ calcAverages[i][j].dist /= (float) iters;
+ }
+ }
+ }else { //find median
+ for (int i = 0; i < calcAverages.size(); i++) { //for each calc
+ for (int j = 0; j < calcAverages[i].size(); j++) { //for each comparison
+ vector<double> dists;
+ for (int thisIter = 0; thisIter < iters; thisIter++) { //for each subsample
+ dists.push_back(calcDistsTotals[thisIter][i][j].dist);
+ }
+ sort(dists.begin(), dists.end());
+ calcAverages[i][j].dist = dists[(iters/2)];
+ }
+ }
+ }
+ //find standard deviation
+ vector< vector<seqDist> > stdDev; stdDev.resize(matrixCalculators.size());
+ for (int i = 0; i < stdDev.size(); i++) { //initialize sums to zero.
+ stdDev[i].resize(calcDistsTotals[0][i].size());
+
+ for (int j = 0; j < stdDev[i].size(); j++) {
+ stdDev[i][j].seq1 = calcDistsTotals[0][i][j].seq1;
+ stdDev[i][j].seq2 = calcDistsTotals[0][i][j].seq2;
+ stdDev[i][j].dist = 0.0;
+ }
+ }
+
+ for (int thisIter = 0; thisIter < iters; thisIter++) { //compute the difference of each dist from the mean, and square the result of each
+ for (int i = 0; i < stdDev.size(); i++) {
+ for (int j = 0; j < stdDev[i].size(); j++) {
+ stdDev[i][j].dist += ((calcDistsTotals[thisIter][i][j].dist - calcAverages[i][j].dist) * (calcDistsTotals[thisIter][i][j].dist - calcAverages[i][j].dist));
+ }
+ }
+ }
+
+ for (int i = 0; i < stdDev.size(); i++) { //finds average.
+ for (int j = 0; j < stdDev[i].size(); j++) {
+ stdDev[i][j].dist /= (float) iters;
+ stdDev[i][j].dist = sqrt(stdDev[i][j].dist);
+ }
+ }
+
+ //print results
+ for (int i = 0; i < calcDists.size(); i++) {
+ vector< vector<double> > matrix; //square matrix to represent the distance
+ matrix.resize(thisLookup.size());
+ for (int k = 0; k < thisLookup.size(); k++) { matrix[k].resize(thisLookup.size(), 0.0); }
+
+ vector< vector<double> > stdmatrix; //square matrix to represent the stdDev
+ stdmatrix.resize(thisLookup.size());
+ for (int k = 0; k < thisLookup.size(); k++) { stdmatrix[k].resize(thisLookup.size(), 0.0); }
+
+
+ for (int j = 0; j < calcAverages[i].size(); j++) {
+ int row = calcAverages[i][j].seq1;
+ int column = calcAverages[i][j].seq2;
+ float dist = calcAverages[i][j].dist;
+ float stdDist = stdDev[i][j].dist;
+
+ matrix[row][column] = dist;
+ matrix[column][row] = dist;
+ stdmatrix[row][column] = stdDist;
+ stdmatrix[column][row] = stdDist;
+ }
+
+ map<string, string> variables;
+ variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(sharedfile));
+ variables["[distance]"] = thisLookup[0]->getLabel();
+ variables["[outputtag]"] = output;
+ variables["[tag2]"] = "ave";
+ variables["[calc]"] = matrixCalculators[i]->getName();
+ string distFileName = getOutputFileName("phylip",variables);
+ outputNames.push_back(distFileName); outputTypes["phylip"].push_back(distFileName);
+ ofstream outAve;
+ m->openOutputFile(distFileName, outAve);
+ outAve.setf(ios::fixed, ios::floatfield); outAve.setf(ios::showpoint);
+
+ printSims(outAve, matrix);
+
+ outAve.close();
+
+ variables["[tag2]"] = "std";
+ distFileName = getOutputFileName("phylip",variables);
+ outputNames.push_back(distFileName); outputTypes["phylip"].push_back(distFileName);
+ ofstream outSTD;
+ m->openOutputFile(distFileName, outSTD);
+ outSTD.setf(ios::fixed, ios::floatfield); outSTD.setf(ios::showpoint);
+
+ printSims(outSTD, stdmatrix);
+
+ outSTD.close();
+
+ }
+ }
+
+ return 0;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "MatrixOutputCommand", "process");
+ exit(1);
+ }
+}
+/**************************************************************************************************/
+int MatrixOutputCommand::driver(vector<SharedRAbundVector*> thisLookup, int start, int end, vector< vector<seqDist> >& calcDists) {
+ try {
+ vector<SharedRAbundVector*> subset;
+
+ for (int k = start; k < end; k++) { // pass cdd each set of groups to compare
+
+ for (int l = 0; l < k; l++) {