]> git.donarmstrong.com Git - mothur.git/blobdiff - unifracweightedcommand.cpp
mods to amova command
[mothur.git] / unifracweightedcommand.cpp
index 9f961a9a6658683a3661290757ebe01048002371..4c387da9332c6ca280a9359105cf8dfa837cd0a9 100644 (file)
@@ -24,11 +24,12 @@ vector<string> UnifracWeightedCommand::getValidParameters(){
 //**********************************************************************************************************************
 UnifracWeightedCommand::UnifracWeightedCommand(){      
        try {
-               //initialize outputTypes
+               abort = true; calledHelp = true; 
                vector<string> tempOutNames;
                outputTypes["weighted"] = tempOutNames;
                outputTypes["wsummary"] = tempOutNames;
                outputTypes["phylip"] = tempOutNames;
+               outputTypes["column"] = tempOutNames;
        }
        catch(exception& e) {
                m->errorOut(e, "UnifracWeightedCommand", "UnifracWeightedCommand");
@@ -63,11 +64,11 @@ vector<string> UnifracWeightedCommand::getRequiredFiles(){
 UnifracWeightedCommand::UnifracWeightedCommand(string option) {
        try {
                globaldata = GlobalData::getInstance();
-               abort = false;
+               abort = false; calledHelp = false;   
                Groups.clear();
                        
                //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
@@ -89,6 +90,7 @@ UnifracWeightedCommand::UnifracWeightedCommand(string option) {
                        outputTypes["weighted"] = tempOutNames;
                        outputTypes["wsummary"] = tempOutNames;
                        outputTypes["phylip"] = tempOutNames;
+                       outputTypes["column"] = tempOutNames;
                        
                        if (globaldata->gTree.size() == 0) {//no trees were read
                                m->mothurOut("You must execute the read.tree command, before you may execute the unifrac.weighted command."); m->mothurOutEndLine(); abort = true;  }
@@ -111,9 +113,13 @@ UnifracWeightedCommand::UnifracWeightedCommand(string option) {
                        itersString = validParameter.validFile(parameters, "iters", false);                     if (itersString == "not found") { itersString = "1000"; }
                        convert(itersString, iters); 
                        
-                       string temp = validParameter.validFile(parameters, "distance", false);                  if (temp == "not found") { temp = "false"; }
-                       phylip = m->isTrue(temp);
-               
+                       string temp = validParameter.validFile(parameters, "distance", false);                  
+                       if (temp == "not found") { phylip = false; outputForm = ""; }
+                       else{
+                               if ((temp == "lt") || (temp == "column") || (temp == "square")) {  phylip = true;  outputForm = temp; }
+                               else { m->mothurOut("Options for distance are: lt, square, or column. Using lt."); m->mothurOutEndLine(); phylip = true; outputForm = "lt"; }
+                       }
+                       
                        temp = validParameter.validFile(parameters, "random", false);                                   if (temp == "not found") { temp = "F"; }
                        random = m->isTrue(temp);
                        
@@ -174,7 +180,7 @@ void UnifracWeightedCommand::help(){
 int UnifracWeightedCommand::execute() {
        try {
        
-               if (abort == true) { return 0; }
+               if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
                
                int start = time(NULL);
                
@@ -211,7 +217,6 @@ int UnifracWeightedCommand::execute() {
                        }
                        
                        if (random) { 
-                               vector<double> sums = weighted->getBranchLengthSums(T[i]); 
                        
                                //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
                                vector< vector<string> > namesOfGroupCombos;
@@ -245,18 +250,18 @@ int UnifracWeightedCommand::execute() {
                                
                                        #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
                                                if(processors == 1){
-                                                       driver(T[i],  namesOfGroupCombos, 0, namesOfGroupCombos.size(), sums, rScores);
+                                                       driver(T[i],  namesOfGroupCombos, 0, namesOfGroupCombos.size(),  rScores);
                                                }else{
-                                                       createProcesses(T[i],  namesOfGroupCombos, sums, rScores);
+                                                       createProcesses(T[i],  namesOfGroupCombos, rScores);
                                                }
                                        #else
-                                               driver(T[i], namesOfGroupCombos, 0, namesOfGroupCombos.size(), sums, rScores);
+                                               driver(T[i], namesOfGroupCombos, 0, namesOfGroupCombos.size(), rScores);
                                        #endif
                                        
                                        if (m->control_pressed) { delete output; outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str());  } return 0; }
                                        
                                        //report progress
-                                       m->mothurOut("Iter: " + toString(j+1)); m->mothurOutEndLine();          
+//                                     m->mothurOut("Iter: " + toString(j+1)); m->mothurOutEndLine();          
                                }
                                lines.clear();
                        
@@ -322,11 +327,10 @@ int UnifracWeightedCommand::execute() {
 }
 /**************************************************************************************************/
 
-int UnifracWeightedCommand::createProcesses(Tree* t, vector< vector<string> > namesOfGroupCombos, vector<double>& sums, vector< vector<double> >& scores) {
+int UnifracWeightedCommand::createProcesses(Tree* t, vector< vector<string> > namesOfGroupCombos, vector< vector<double> >& scores) {
        try {
 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
                int process = 1;
-               int num = 0;
                vector<int> processIDS;
                
                EstOutput results;
@@ -339,7 +343,7 @@ int UnifracWeightedCommand::createProcesses(Tree* t, vector< vector<string> > na
                                processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
                                process++;
                        }else if (pid == 0){
-                               driver(t, namesOfGroupCombos, lines[process].start, lines[process].num, sums, scores);
+                               driver(t, namesOfGroupCombos, lines[process].start, lines[process].num, scores);
                        
                                //pass numSeqs to parent
                                ofstream out;
@@ -349,10 +353,14 @@ int UnifracWeightedCommand::createProcesses(Tree* t, vector< vector<string> > na
                                out.close();
                                
                                exit(0);
-                       }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
+                       }else { 
+                               m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); 
+                               for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
+                               exit(0);
+                       }
                }
                
-               driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, sums, scores);
+               driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, scores);
                
                //force parent to wait until all the processes are done
                for (int i=0;i<(processors-1);i++) { 
@@ -383,7 +391,7 @@ int UnifracWeightedCommand::createProcesses(Tree* t, vector< vector<string> > na
 }
 
 /**************************************************************************************************/
-int UnifracWeightedCommand::driver(Tree* t, vector< vector<string> > namesOfGroupCombos, int start, int num, vector<double>& sums, vector< vector<double> >& scores) { 
+int UnifracWeightedCommand::driver(Tree* t, vector< vector<string> > namesOfGroupCombos, int start, int num, vector< vector<double> >& scores) { 
  try {
                Tree* randT = new Tree();
 
@@ -404,7 +412,7 @@ int UnifracWeightedCommand::driver(Tree* t, vector< vector<string> > namesOfGrou
                        if (m->control_pressed) { delete randT;  return 0;  }
 
                        //get wscore of random tree
-                       EstOutput randomData = weighted->getValues(randT, groupA, groupB, sums);
+                       EstOutput randomData = weighted->getValues(randT, groupA, groupB);
                
                        if (m->control_pressed) { delete randT;  return 0;  }
                                                                                
@@ -493,15 +501,23 @@ void UnifracWeightedCommand::createPhylipFile() {
                //for each tree
                for (int i = 0; i < T.size(); i++) { 
                
-                       string phylipFileName = outputDir + m->getSimpleName(globaldata->getTreeFile())  + toString(i+1) + ".weighted.dist";
-                       outputNames.push_back(phylipFileName);
-                       outputTypes["phylip"].push_back(phylipFileName);
+                       string phylipFileName;
+                       if ((outputForm == "lt") || (outputForm == "square")) {
+                               phylipFileName = outputDir + m->getSimpleName(globaldata->getTreeFile())  + toString(i+1) + ".weighted.phylip.dist";
+                               outputNames.push_back(phylipFileName); outputTypes["phylip"].push_back(phylipFileName); 
+                       }else { //column
+                               phylipFileName = outputDir + m->getSimpleName(globaldata->getTreeFile())  + toString(i+1) + ".weighted.column.dist";
+                               outputNames.push_back(phylipFileName); outputTypes["column"].push_back(phylipFileName); 
+                       }
+                       
                        ofstream out;
                        m->openOutputFile(phylipFileName, out);
                        
-                       //output numSeqs
-                       out << globaldata->Groups.size() << endl;
-                       
+                       if ((outputForm == "lt") || (outputForm == "square")) {
+                               //output numSeqs
+                               out << globaldata->Groups.size() << endl;
+                       }
+
                        //make matrix with scores in it
                        vector< vector<float> > dists;  dists.resize(globaldata->Groups.size());
                        for (int i = 0; i < globaldata->Groups.size(); i++) {
@@ -510,7 +526,7 @@ void UnifracWeightedCommand::createPhylipFile() {
                        
                        //flip it so you can print it
                        for (int r=0; r<globaldata->Groups.size(); r++) { 
-                               for (int l = r+1; l < globaldata->Groups.size(); l++) {
+                               for (int l = 0; l < r; l++) {
                                        dists[r][l] = utreeScores[count];
                                        dists[l][r] = utreeScores[count];
                                        count++;
@@ -524,11 +540,30 @@ void UnifracWeightedCommand::createPhylipFile() {
                                if (name.length() < 10) { //pad with spaces to make compatible
                                        while (name.length() < 10) {  name += " ";  }
                                }
-                               out << name << '\t';
                                
-                               //output distances
-                               for (int l = 0; l < r; l++) {   out  << dists[r][l] << '\t';  }
-                               out << endl;
+                               if (outputForm == "lt") {
+                                       out << name << '\t';
+                                       
+                                       //output distances
+                                       for (int l = 0; l < r; l++) {   out  << dists[r][l] << '\t';  }
+                                       out << endl;
+                               }else if (outputForm == "square") {
+                                       out << name << '\t';
+                                       
+                                       //output distances
+                                       for (int l = 0; l < globaldata->Groups.size(); l++) {   out  << dists[r][l] << '\t';  }
+                                       out << endl;
+                               }else{
+                                       //output distances
+                                       for (int l = 0; l < r; l++) {   
+                                               string otherName = globaldata->Groups[l];
+                                               if (otherName.length() < 10) { //pad with spaces to make compatible
+                                                       while (otherName.length() < 10) {  otherName += " ";  }
+                                               }
+                                               
+                                               out  << name << '\t' << otherName << dists[r][l] << endl;  
+                                       }
+                               }
                        }
                        out.close();
                }