changed reading of name file to use buffered reads. note the splitAtWhiteSpace functi...

[mothur.git] / hclustercommand.cpp

diff --git a/hclustercommand.cpp b/hclustercommand.cpp
index 582fc45f3f4744caf50705e6a10482509fbd4afe..47136938a19e6c5d4b8dd2b8110849249ac50b06 100644 (file)
--- a/hclustercommand.cpp
+++ b/hclustercommand.cpp
@@ -17,8 +17,8 @@ vector<string> HClusterCommand::setParameters(){
                CommandParameter pcolumn("column", "InputTypes", "", "", "PhylipColumn", "PhylipColumn", "ColumnName",false,false); parameters.push_back(pcolumn);
                CommandParameter pcutoff("cutoff", "Number", "", "10", "", "", "",false,false); parameters.push_back(pcutoff);
                CommandParameter pprecision("precision", "Number", "", "100", "", "", "",false,false); parameters.push_back(pprecision);
-               CommandParameter pmethod("method", "Multiple", "furthest-nearest-average-weighted", "furthest", "", "", "",false,false); parameters.push_back(pmethod);
-               CommandParameter phard("hard", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(phard);
+               CommandParameter pmethod("method", "Multiple", "furthest-nearest-average-weighted", "average", "", "", "",false,false); parameters.push_back(pmethod);
+               CommandParameter phard("hard", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(phard);
                CommandParameter psorted("sorted", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(psorted);
                CommandParameter pshowabund("showabund", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pshowabund);
                CommandParameter ptiming("timing", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(ptiming);          
@@ -43,7 +43,7 @@ string HClusterCommand::getHelpString(){
                helpString += "The name parameter allows you to enter your name file and is required if your distance file is in column format. \n";
                helpString += "The hcluster command should be in the following format: \n";
                helpString += "hcluster(column=youDistanceFile, name=yourNameFile, method=yourMethod, cutoff=yourCutoff, precision=yourPrecision) \n";
-               helpString += "The acceptable hcluster methods are furthest, nearest, weighted and average.\n\n";       
+               helpString += "The acceptable hcluster methods are furthest, nearest, weighted and average.\n"; 
                return helpString;
        }
        catch(exception& e) {
@@ -74,11 +74,10 @@ HClusterCommand::HClusterCommand(string option)  {
                
                //allow user to run help
                if(option == "help") { help(); abort = true; calledHelp = true; }
+               else if(option == "citation") { citation(); abort = true; calledHelp = true;}
                
                else {
-                       //valid paramters for this command
-                       string Array[] =  {"cutoff","hard","precision","method","phylip","column","name","sorted","showabund","timing","outputdir","inputdir"};
-                       vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+                       vector<string> myArray = setParameters();
                        
                        OptionParser parser(option);
                        map<string,string> parameters = parser.getParameters();
@@ -136,16 +135,17 @@ HClusterCommand::HClusterCommand(string option)  {
                        phylipfile = validParameter.validFile(parameters, "phylip", true);
                        if (phylipfile == "not open") { abort = true; }
                        else if (phylipfile == "not found") { phylipfile = ""; }        
-                       else {  distfile = phylipfile;  format = "phylip";      }
+                       else {  distfile = phylipfile;  format = "phylip";      m->setPhylipFile(phylipfile); }
                        
                        columnfile = validParameter.validFile(parameters, "column", true);
                        if (columnfile == "not open") { abort = true; } 
                        else if (columnfile == "not found") { columnfile = ""; }
-                       else {  distfile = columnfile; format = "column";       }
+                       else {  distfile = columnfile; format = "column";       m->setColumnFile(columnfile); }
                        
                        namefile = validParameter.validFile(parameters, "name", true);
                        if (namefile == "not open") { abort = true; }   
                        else if (namefile == "not found") { namefile = ""; }
+                       else { m->setNameFile(namefile); }
                        
                        if ((phylipfile == "") && (columnfile == "")) { 
                                //is there are current file available for either of these?
@@ -182,18 +182,18 @@ HClusterCommand::HClusterCommand(string option)  {
                        if (temp == "not found") { temp = "100"; }
                        //saves precision legnth for formatting below
                        length = temp.length();
-                       convert(temp, precision); 
+                       m->mothurConvert(temp, precision); 
                        
-                       temp = validParameter.validFile(parameters, "hard", false);                     if (temp == "not found") { temp = "F"; }
+                       temp = validParameter.validFile(parameters, "hard", false);                     if (temp == "not found") { temp = "T"; }
                        hard = m->isTrue(temp);
                        
                        temp = validParameter.validFile(parameters, "cutoff", false);
                        if (temp == "not found") { temp = "10"; }
-                       convert(temp, cutoff); 
+                       m->mothurConvert(temp, cutoff); 
                        cutoff += (5 / (precision * 10.0)); 
                        
                        method = validParameter.validFile(parameters, "method", false);
-                       if (method == "not found") { method = "furthest"; }
+                       if (method == "not found") { method = "average"; }
                        
                        if ((method == "furthest") || (method == "nearest") || (method == "average") || (method == "weighted")) { }
                        else { m->mothurOut("Not a valid clustering method.  Valid clustering algorithms are furthest, nearest, average or weighted."); m->mothurOutEndLine(); abort = true; }
@@ -260,7 +260,7 @@ int HClusterCommand::execute(){
                                sabundFile.close();
                                rabundFile.close();
                                listFile.close();
-                               for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  } outputTypes.clear();
+                               for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  } outputTypes.clear();
                                return 0;  
                        }
                        
@@ -276,7 +276,7 @@ int HClusterCommand::execute(){
                        sabundFile.close();
                        rabundFile.close();
                        listFile.close();
-                       for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  } outputTypes.clear();
+                       for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  } outputTypes.clear();
                        return 0;  
                }
 
@@ -306,35 +306,41 @@ int HClusterCommand::execute(){
                                sabundFile.close();
                                rabundFile.close();
                                listFile.close();
-                               for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  } outputTypes.clear();
+                               for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  } outputTypes.clear();
                                return 0;  
                }
 
+               float saveCutoff = cutoff;
                
                while (seqs.size() != 0){
                
                        seqs = cluster->getSeqs();
                        
+                       //to account for cutoff change in average neighbor
+                       if (seqs.size() != 0) {
+                               if (seqs[0].dist > cutoff) { break; }
+                       }
+                       
                        if (m->control_pressed) {  
                                delete cluster;
                                sabundFile.close();
                                rabundFile.close();
                                listFile.close();
-                               for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  } outputTypes.clear();
+                               for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  } outputTypes.clear();
                                return 0;  
                        }
 
                        for (int i = 0; i < seqs.size(); i++) {  //-1 means skip me
                                
                                if (seqs[i].seq1 != seqs[i].seq2) {
-                                       cluster->update(seqs[i].seq1, seqs[i].seq2, seqs[i].dist);
+                                       cutoff = cluster->update(seqs[i].seq1, seqs[i].seq2, seqs[i].dist);
                                        
                                        if (m->control_pressed) {  
                                                delete cluster;
                                                sabundFile.close();
                                                rabundFile.close();
                                                listFile.close();
-                                               for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  } outputTypes.clear();
+                                               for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  } outputTypes.clear();
                                                return 0;  
                                        }
 
@@ -367,7 +373,7 @@ int HClusterCommand::execute(){
                        sabundFile.close();
                        rabundFile.close();
                        listFile.close();
-                       for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  } outputTypes.clear();
+                       for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  } outputTypes.clear();
                        return 0;  
                }
                                        
@@ -384,10 +390,18 @@ int HClusterCommand::execute(){
                delete cluster;
                
                if (m->control_pressed) {  
-                       for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  } outputTypes.clear();
+                       for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  } outputTypes.clear();
                        return 0;  
                }
                
+               
+               if (saveCutoff != cutoff) { 
+                       if (hard)       {  saveCutoff = m->ceilDist(saveCutoff, precision);     }
+                       else            {       saveCutoff = m->roundDist(saveCutoff, precision);  }
+                       
+                       m->mothurOut("changed cutoff to " + toString(cutoff)); m->mothurOutEndLine(); 
+               }
+               
                //set list file as new current listfile
                string current = "";
                itTypes = outputTypes.find("list");