some changes to trim.flows

[mothur.git] / nmdscommand.cpp

diff --git a/nmdscommand.cpp b/nmdscommand.cpp
index 4adb57c6147c544ab9284996d85a8bf4db9e55e2..aa1cd660b9af25b9dfdb099741b06fb6769433c1 100644 (file)
--- a/nmdscommand.cpp
+++ b/nmdscommand.cpp
@@ -25,8 +25,7 @@ vector<string> NMDSCommand::getValidParameters(){
 //**********************************************************************************************************************
 NMDSCommand::NMDSCommand(){    
        try {
-               abort = true;
-               //initialize outputTypes
+               abort = true; calledHelp = true; 
                vector<string> tempOutNames;
                outputTypes["nmds"] = tempOutNames;
                outputTypes["stress"] = tempOutNames;
@@ -64,10 +63,10 @@ vector<string> NMDSCommand::getRequiredFiles(){
 
 NMDSCommand::NMDSCommand(string option)  {
        try {
-               abort = false;
+               abort = false; calledHelp = false;   
                
                //allow user to run help
-               if(option == "help") { help(); abort = true; }
+               if(option == "help") { help(); abort = true; calledHelp = true; }
                
                else {
                        //valid paramters for this command
@@ -127,7 +126,7 @@ NMDSCommand::NMDSCommand(string option)  {
                                outputDir += m->hasPath(phylipfile); //if user entered a file with a path then preserve it      
                        }
                        
-                       string temp = validParameter.validFile(parameters, "mindim", false);    if (temp == "not found") {      temp = "1";     }
+                       string temp = validParameter.validFile(parameters, "mindim", false);    if (temp == "not found") {      temp = "2";     }
                        convert(temp, mindim);
                        
                        temp = validParameter.validFile(parameters, "maxiters", false); if (temp == "not found") {      temp = "500";   }
@@ -159,8 +158,8 @@ void NMDSCommand::help(){
                m->mothurOut("The nmds command parameters are phylip, axes, mindim, maxdim, maxiters, iters and epsilon."); m->mothurOutEndLine();
                m->mothurOut("The phylip parameter allows you to enter your distance file."); m->mothurOutEndLine();
                m->mothurOut("The axes parameter allows you to enter a file containing a starting configuration."); m->mothurOutEndLine();
-               m->mothurOut("The maxdim parameter allows you to select how maximum dimensions to use. Default=2"); m->mothurOutEndLine();
-               m->mothurOut("The mindim parameter allows you to select how minimum dimensions to use. Default=1"); m->mothurOutEndLine();
+               m->mothurOut("The maxdim parameter allows you to select the maximum dimensions to use. Default=2"); m->mothurOutEndLine();
+               m->mothurOut("The mindim parameter allows you to select the minimum dimensions to use. Default=2"); m->mothurOutEndLine();
                m->mothurOut("The maxiters parameter allows you to select the maximum number of iters to try with each random configuration. Default=500"); m->mothurOutEndLine();
                m->mothurOut("The iters parameter allows you to select the number of random configuration to try. Default=10"); m->mothurOutEndLine();
                m->mothurOut("The epsilon parameter allows you to select set an acceptable stopping point. Default=1e-12."); m->mothurOutEndLine();
@@ -178,7 +177,7 @@ NMDSCommand::~NMDSCommand(){}
 int NMDSCommand::execute(){
        try {
                
-               if (abort == true) { return 0; }
+               if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
                
                cout.setf(ios::fixed, ios::floatfield);
                cout.setf(ios::showpoint);
@@ -196,7 +195,7 @@ int NMDSCommand::execute(){
                if (axesfile != "") {  axes = readAxes(names);          }
                
                string outputFileName = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "nmds.iters";
-               string stressFileName = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "stress.nmds";
+               string stressFileName = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "nmds.stress";
                outputNames.push_back(outputFileName); outputTypes["iters"].push_back(outputFileName);
                outputNames.push_back(stressFileName); outputTypes["stress"].push_back(stressFileName);
                
@@ -209,10 +208,12 @@ int NMDSCommand::execute(){
                out.setf(ios::fixed, ios::floatfield);
                out.setf(ios::showpoint);
                
-               out2 << "Dimension\tIter\tStress\tCorr" << endl;
+               out2 << "Dimension\tIter\tStress\tRsq" << endl;
                
                double bestStress = 10000000;
+               double bestR2 = 10000000;
                vector< vector<double> > bestConfig;
+               int bestDim = 0;
                
                for (int i = mindim; i <= maxdim; i++) {
                        m->mothurOut("Processing Dimension: " + toString(i)); m->mothurOutEndLine();
@@ -236,21 +237,23 @@ int NMDSCommand::execute(){
                                if (m->control_pressed) { out.close(); out2.close(); for (int k = 0; k < outputNames.size(); k++) {     remove(outputNames[k].c_str()); } return 0; }
                                
                                //calc correlation between original distances and euclidean distances from this config
-                               double corr = linearCalc.calcPearson(newEuclid, matrix);
-                               corr *= corr;
+                               double rsquared = linearCalc.calcPearson(newEuclid, matrix);
+                               rsquared *= rsquared;
                                if (m->control_pressed) { out.close(); out2.close(); for (int k = 0; k < outputNames.size(); k++) {     remove(outputNames[k].c_str()); } return 0; }
                                
                                //output results
                                out << "Config" << (j+1) << '\t';
                                for (int k = 0; k < i; k++) { out << "axis" << (k+1) << '\t'; }
                                out << endl;
-                               out2 << i << '\t' << (j+1) << '\t' << stress << '\t' << corr << endl;
+                               out2 << i << '\t' << (j+1) << '\t' << stress << '\t' << rsquared << endl;
                                
                                output(endConfig, names, out);
                                
                                //save best
                                if (stress < bestStress) {
+                                       bestDim = i;
                                        bestStress = stress;
+                                       bestR2 = rsquared;
                                        bestConfig = endConfig;
                                }
                                
@@ -264,6 +267,10 @@ int NMDSCommand::execute(){
                string BestFileName = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "nmds.axes";
                outputNames.push_back(BestFileName); outputTypes["nmds"].push_back(BestFileName);
                
+               m->mothurOut("\nNumber of dimensions:\t" + toString(bestDim) + "\n");
+               m->mothurOut("Lowest stress :\t" + toString(bestStress) + "\n");
+               m->mothurOut("R-squared for configuration:\t" + toString(bestR2) + "\n");
+               
                ofstream outBest;
                m->openOutputFile(BestFileName, outBest);
                outBest.setf(ios::fixed, ios::floatfield);