parsimonycommand.cpp

   1 /*
   2  *  parsimonycommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 1/26/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "parsimonycommand.h"
  11
  12 //**********************************************************************************************************************
  13 vector<string> ParsimonyCommand::setParameters(){
  14         try {
  15                 CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(ptree);
  16                 CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pgroup);
  17                 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
  18                 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
  19                 CommandParameter prandom("random", "String", "", "", "", "", "",false,false); parameters.push_back(prandom);
  20                 CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters);
  21                 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
  22                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  23                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  24
  25                 vector<string> myArray;
  26                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  27                 return myArray;
  28         }
  29         catch(exception& e) {
  30                 m->errorOut(e, "ParsimonyCommand", "setParameters");
  31                 exit(1);
  32         }
  33 }
  34 //**********************************************************************************************************************
  35 string ParsimonyCommand::getHelpString(){
  36         try {
  37                 string helpString = "";
  38                 helpString += "The parsimony command parameters are tree, group, name, random, groups, processors and iters.  tree parameter is required unless you have valid current tree file or are using random.\n";
  39                 helpString += "The groups parameter allows you to specify which of the groups in your groupfile you would like analyzed.  You must enter at least 1 valid group.\n";
  40                 helpString += "The group names are separated by dashes.  The iters parameter allows you to specify how many random trees you would like compared to your tree.\n";
  41                 helpString += "The parsimony command should be in the following format: parsimony(random=yourOutputFilename, groups=yourGroups, iters=yourIters).\n";
  42                 helpString += "The processors parameter allows you to specify the number of processors to use. The default is 1.\n";
  43                 helpString += "Example parsimony(random=out, iters=500).\n";
  44                 helpString += "The default value for random is "" (meaning you want to use the trees in your inputfile, randomtree=out means you just want the random distribution of trees outputted to out.rd_parsimony),\n";
  45                 helpString += "and iters is 1000.  The parsimony command output two files: .parsimony and .psummary their descriptions are in the manual.\n";
  46                 helpString += "Note: No spaces between parameter labels (i.e. random), '=' and parameters (i.e.yourOutputFilename).\n";
  47                 return helpString;
  48         }
  49         catch(exception& e) {
  50                 m->errorOut(e, "ParsimonyCommand", "getHelpString");
  51                 exit(1);
  52         }
  53 }
  54
  55 //**********************************************************************************************************************
  56 ParsimonyCommand::ParsimonyCommand(){
  57         try {
  58                 abort = true; calledHelp = true;
  59                 setParameters();
  60                 vector<string> tempOutNames;
  61                 outputTypes["parsimony"] = tempOutNames;
  62                 outputTypes["psummary"] = tempOutNames;
  63         }
  64         catch(exception& e) {
  65                 m->errorOut(e, "ParsimonyCommand", "ParsimonyCommand");
  66                 exit(1);
  67         }
  68 }
  69 /***********************************************************/
  70 ParsimonyCommand::ParsimonyCommand(string option)  {
  71         try {
  72                 abort = false; calledHelp = false;
  73                 Groups.clear();
  74
  75                 //allow user to run help
  76                 if(option == "help") { help(); abort = true; calledHelp = true; }
  77                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  78
  79                 else {
  80                         vector<string> myArray = setParameters();
  81
  82                         OptionParser parser(option);
  83                         map<string, string> parameters = parser.getParameters();
  84                         map<string,string>::iterator it;
  85
  86                         ValidParameters validParameter;
  87
  88                         //check to make sure all parameters are valid for command
  89                         for (it = parameters.begin(); it != parameters.end(); it++) {
  90                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  91                         }
  92
  93                         //initialize outputTypes
  94                         vector<string> tempOutNames;
  95                         outputTypes["parsimony"] = tempOutNames;
  96                         outputTypes["psummary"] = tempOutNames;
  97
  98                         //if the user changes the input directory command factory will send this info to us in the output parameter
  99                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 100                         if (inputDir == "not found"){   inputDir = "";          }
 101                         else {
 102                                 string path;
 103                                 it = parameters.find("tree");
 104                                 //user has given a template file
 105                                 if(it != parameters.end()){
 106                                         path = m->hasPath(it->second);
 107                                         //if the user has not given a path then, add inputdir. else leave path alone.
 108                                         if (path == "") {       parameters["tree"] = inputDir + it->second;             }
 109                                 }
 110
 111                                 it = parameters.find("group");
 112                                 //user has given a template file
 113                                 if(it != parameters.end()){
 114                                         path = m->hasPath(it->second);
 115                                         //if the user has not given a path then, add inputdir. else leave path alone.
 116                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
 117                                 }
 118
 119                                 it = parameters.find("name");
 120                                 //user has given a template file
 121                                 if(it != parameters.end()){
 122                                         path = m->hasPath(it->second);
 123                                         //if the user has not given a path then, add inputdir. else leave path alone.
 124                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
 125                                 }
 126                         }
 127
 128                         m->runParse = true;
 129                         m->clearGroups();
 130                         m->clearAllGroups();
 131                         m->Treenames.clear();
 132                         m->names.clear();
 133
 134                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = ""; }
 135
 136                         randomtree = validParameter.validFile(parameters, "random", false);             if (randomtree == "not found") { randomtree = ""; }
 137
 138                         //are you trying to use parsimony without reading a tree or saying you want random distribution
 139                         if (randomtree == "")  {
 140                                 //check for required parameters
 141                                 treefile = validParameter.validFile(parameters, "tree", true);
 142                                 if (treefile == "not open") { treefile = ""; abort = true; }
 143                                 else if (treefile == "not found") {                             //if there is a current design file, use it
 144                                         treefile = m->getTreeFile();
 145                                         if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
 146                                         else {  m->mothurOut("You have no current tree file and the tree parameter is required."); m->mothurOutEndLine(); abort = true; }
 147                                 }else { m->setTreeFile(treefile); }
 148
 149                                 //check for required parameters
 150                                 groupfile = validParameter.validFile(parameters, "group", true);
 151                                 if (groupfile == "not open") { abort = true; }
 152                                 else if (groupfile == "not found") { groupfile = ""; }
 153                                 else { m->setGroupFile(groupfile); }
 154
 155                                 namefile = validParameter.validFile(parameters, "name", true);
 156                                 if (namefile == "not open") { namefile = ""; abort = true; }
 157                                 else if (namefile == "not found") { namefile = ""; }
 158                                 else { m->setNameFile(namefile); }
 159                         }
 160
 161                         //if the user changes the output directory command factory will send this info to us in the output parameter
 162                         string outputDir = validParameter.validFile(parameters, "outputdir", false);            if (outputDir == "not found"){  outputDir = ""; if (randomtree == "")  { outputDir += m->hasPath(treefile); } }
 163
 164                         //check for optional parameter and set defaults
 165                         // ...at some point should added some additional type checking...
 166                         groups = validParameter.validFile(parameters, "groups", false);
 167                         if (groups == "not found") { groups = ""; m->clearGroups(); }
 168                         else {
 169                                 m->splitAtDash(groups, Groups);
 170                                 m->setGroups(Groups);
 171                         }
 172
 173                         itersString = validParameter.validFile(parameters, "iters", false);                     if (itersString == "not found") { itersString = "1000"; }
 174                         m->mothurConvert(itersString, iters);
 175
 176                         string temp = validParameter.validFile(parameters, "processors", false);        if (temp == "not found"){       temp = m->getProcessors();      }
 177                         m->setProcessors(temp);
 178                         m->mothurConvert(temp, processors);
 179
 180                         if (namefile == "") {
 181                                 vector<string> files; files.push_back(treefile);
 182                                 parser.getNameFile(files);
 183                         }
 184
 185                 }
 186
 187         }
 188         catch(exception& e) {
 189                 m->errorOut(e, "ParsimonyCommand", "ParsimonyCommand");
 190                 exit(1);
 191         }
 192 }
 193 /***********************************************************/
 194 int ParsimonyCommand::execute() {
 195         try {
 196
 197                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 198
 199
 200                 //randomtree will tell us if user had their own treefile or if they just want the random distribution
 201                 //user has entered their own tree
 202                 if (randomtree == "") {
 203
 204                         m->setTreeFile(treefile);
 205
 206                         if (groupfile != "") {
 207                                 //read in group map info.
 208                                 tmap = new TreeMap(groupfile);
 209                                 tmap->readMap();
 210                         }else{ //fake out by putting everyone in one group
 211                                 Tree* tree = new Tree(treefile); delete tree;  //extracts names from tree to make faked out groupmap
 212                                 tmap = new TreeMap();
 213
 214                                 for (int i = 0; i < m->Treenames.size(); i++) { tmap->addSeq(m->Treenames[i], "Group1"); }
 215                         }
 216
 217                         if (namefile != "") { readNamesFile(); }
 218
 219                         read = new ReadNewickTree(treefile);
 220                         int readOk = read->read(tmap);
 221
 222                         if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete tmap; delete read; return 0; }
 223
 224                         read->AssembleTrees();
 225                         T = read->getTrees();
 226                         delete read;
 227
 228                         //make sure all files match
 229                         //if you provide a namefile we will use the numNames in the namefile as long as the number of unique match the tree names size.
 230                         int numNamesInTree;
 231                         if (namefile != "")  {
 232                                 if (numUniquesInName == m->Treenames.size()) {  numNamesInTree = nameMap.size();  }
 233                                 else {   numNamesInTree = m->Treenames.size();  }
 234                         }else {  numNamesInTree = m->Treenames.size();  }
 235
 236
 237                         //output any names that are in group file but not in tree
 238                         if (numNamesInTree < tmap->getNumSeqs()) {
 239                                 for (int i = 0; i < tmap->namesOfSeqs.size(); i++) {
 240                                         //is that name in the tree?
 241                                         int count = 0;
 242                                         for (int j = 0; j < m->Treenames.size(); j++) {
 243                                                 if (tmap->namesOfSeqs[i] == m->Treenames[j]) { break; } //found it
 244                                                 count++;
 245                                         }
 246
 247                                         if (m->control_pressed) {
 248                                                 delete tmap; for (int i = 0; i < T.size(); i++) { delete T[i]; }
 249                                                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]); } outputTypes.clear();
 250                                                 m->clearGroups();
 251                                                 return 0;
 252                                         }
 253
 254                                         //then you did not find it so report it
 255                                         if (count == m->Treenames.size()) {
 256                                                 //if it is in your namefile then don't remove
 257                                                 map<string, string>::iterator it = nameMap.find(tmap->namesOfSeqs[i]);
 258
 259                                                 if (it == nameMap.end()) {
 260                                                         m->mothurOut(tmap->namesOfSeqs[i] + " is in your groupfile and not in your tree. It will be disregarded."); m->mothurOutEndLine();
 261                                                         tmap->removeSeq(tmap->namesOfSeqs[i]);
 262                                                         i--; //need this because removeSeq removes name from namesOfSeqs
 263                                                 }
 264                                         }
 265                                 }
 266                         }
 267
 268                         if(outputDir == "") { outputDir += m->hasPath(treefile); }
 269                         output = new ColumnFile(outputDir + m->getSimpleName(treefile)  +  ".parsimony", itersString);
 270                         outputNames.push_back(outputDir + m->getSimpleName(treefile)  +  ".parsimony");
 271                         outputTypes["parsimony"].push_back(outputDir + m->getSimpleName(treefile)  +  ".parsimony");
 272
 273                         sumFile = outputDir + m->getSimpleName(treefile) + ".psummary";
 274                         m->openOutputFile(sumFile, outSum);
 275                         outputNames.push_back(sumFile);
 276                         outputTypes["psummary"].push_back(sumFile);
 277                 }else { //user wants random distribution
 278                         getUserInput();
 279
 280                         if(outputDir == "") { outputDir += m->hasPath(randomtree); }
 281                         output = new ColumnFile(outputDir+ m->getSimpleName(randomtree), itersString);
 282                         outputNames.push_back(outputDir+ m->getSimpleName(randomtree));
 283                         outputTypes["parsimony"].push_back(outputDir+ m->getSimpleName(randomtree));
 284                 }
 285
 286                 //set users groups to analyze
 287                 util = new SharedUtil();
 288                 vector<string> mGroups = m->getGroups();
 289                 vector<string> tGroups = tmap->getNamesOfGroups();
 290                 util->setGroups(mGroups, tGroups, allGroups, numGroups, "parsimony");   //sets the groups the user wants to analyze
 291                 util->getCombos(groupComb, mGroups, numComp);
 292                 m->setGroups(mGroups);
 293                 delete util;
 294
 295                 if (numGroups == 1) { numComp++; groupComb.push_back(allGroups); }
 296
 297                 pars = new Parsimony(tmap);
 298                 counter = 0;
 299
 300                 Progress* reading;
 301                 reading = new Progress("Comparing to random:", iters);
 302
 303                 if (m->control_pressed) {
 304                         delete reading; delete pars; delete output;
 305                         delete tmap; for (int i = 0; i < T.size(); i++) { delete T[i]; }
 306                         if (randomtree == "") {  outSum.close();  }
 307                         for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]); } outputTypes.clear();
 308                         m->clearGroups();
 309                         return 0;
 310                 }
 311
 312
 313                 //get pscore for users tree
 314                 userData.resize(numComp,0);  //data = AB, AC, BC, ABC.
 315                 randomData.resize(numComp,0);  //data = AB, AC, BC, ABC.
 316                 rscoreFreq.resize(numComp);
 317                 uscoreFreq.resize(numComp);
 318                 rCumul.resize(numComp);
 319                 uCumul.resize(numComp);
 320                 userTreeScores.resize(numComp);
 321                 UScoreSig.resize(numComp);
 322
 323                 if (randomtree == "") {
 324                         //get pscores for users trees
 325                         for (int i = 0; i < T.size(); i++) {
 326                                 userData = pars->getValues(T[i], processors, outputDir);  //data = AB, AC, BC, ABC.
 327
 328                                 if (m->control_pressed) {
 329                                         delete reading; delete pars; delete output;
 330                                         delete tmap; for (int i = 0; i < T.size(); i++) { delete T[i]; }
 331                                         if (randomtree == "") {  outSum.close();  }
 332                                         for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]); } outputTypes.clear();
 333                                         m->clearGroups();
 334                                         return 0;
 335                                 }
 336
 337
 338                                 //output scores for each combination
 339                                 for(int k = 0; k < numComp; k++) {
 340
 341                                         //update uscoreFreq
 342                                         map<int,double>::iterator it = uscoreFreq[k].find(userData[k]);
 343                                         if (it == uscoreFreq[k].end()) {//new score
 344                                                 uscoreFreq[k][userData[k]] = 1;
 345                                         }else{ uscoreFreq[k][userData[k]]++; }
 346
 347                                         //add users score to valid scores
 348                                         validScores[userData[k]] = userData[k];
 349
 350                                         //save score for summary file
 351                                         userTreeScores[k].push_back(userData[k]);
 352                                 }
 353                         }
 354
 355                         //get pscores for random trees
 356                         for (int j = 0; j < iters; j++) {
 357
 358                                 //create new tree with same num nodes and leaves as users
 359                                 randT = new Tree(tmap);
 360
 361                                 //create random relationships between nodes
 362                                 randT->assembleRandomTree();
 363
 364                                 //get pscore of random tree
 365                                 randomData = pars->getValues(randT, processors, outputDir);
 366
 367                                 if (m->control_pressed) {
 368                                         delete reading; delete pars; delete output; delete randT;
 369                                         if (randomtree == "") {  outSum.close();  }
 370                                         for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]); } outputTypes.clear();
 371                                         delete tmap; for (int i = 0; i < T.size(); i++) { delete T[i]; }
 372                                         m->clearGroups();
 373                                         return 0;
 374                                 }
 375
 376                                 for(int r = 0; r < numComp; r++) {
 377                                         //add trees pscore to map of scores
 378                                         map<int,double>::iterator it = rscoreFreq[r].find(randomData[r]);
 379                                         if (it != rscoreFreq[r].end()) {//already have that score
 380                                                 rscoreFreq[r][randomData[r]]++;
 381                                         }else{//first time we have seen this score
 382                                                 rscoreFreq[r][randomData[r]] = 1;
 383                                         }
 384
 385                                         //add randoms score to validscores
 386                                         validScores[randomData[r]] = randomData[r];
 387                                 }
 388
 389                                 //update progress bar
 390                                 reading->update(j);
 391
 392                                 delete randT;
 393                         }
 394
 395                 }else {
 396                         //get pscores for random trees
 397                         for (int j = 0; j < iters; j++) {
 398
 399                                 //create new tree with same num nodes and leaves as users
 400                                 randT = new Tree(tmap);
 401                                 //create random relationships between nodes
 402
 403                                 randT->assembleRandomTree();
 404
 405                                 if (m->control_pressed) {
 406                                         delete reading; delete pars; delete output; delete randT;
 407                                         delete tmap;
 408                                         for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]); } outputTypes.clear();
 409                                         m->clearGroups();
 410                                         return 0;
 411                                 }
 412
 413
 414                                 //get pscore of random tree
 415                                 randomData = pars->getValues(randT, processors, outputDir);
 416
 417                                 if (m->control_pressed) {
 418                                         delete reading; delete pars;  delete output; delete randT;
 419                                         delete tmap;
 420                                         for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]); } outputTypes.clear();
 421                                         m->clearGroups();
 422                                         return 0;
 423                                 }
 424
 425                                 for(int r = 0; r < numComp; r++) {
 426                                         //add trees pscore to map of scores
 427                                         map<int,double>::iterator it = rscoreFreq[r].find(randomData[r]);
 428                                         if (it != rscoreFreq[r].end()) {//already have that score
 429                                                 rscoreFreq[r][randomData[r]]++;
 430                                         }else{//first time we have seen this score
 431                                                 rscoreFreq[r][randomData[r]] = 1;
 432                                         }
 433
 434                                         //add randoms score to validscores
 435                                         validScores[randomData[r]] = randomData[r];
 436                                 }
 437
 438                                 //update progress bar
 439                                 reading->update(j);
 440
 441                                 delete randT;
 442                         }
 443                 }
 444
 445                 for(int a = 0; a < numComp; a++) {
 446                         float rcumul = 0.0000;
 447                         float ucumul = 0.0000;
 448                         //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
 449                         for (map<int,double>::iterator it = validScores.begin(); it != validScores.end(); it++) {
 450                                 if (randomtree == "") {
 451                                         map<int,double>::iterator it2 = uscoreFreq[a].find(it->first);
 452                                         //user data has that score
 453                                         if (it2 != uscoreFreq[a].end()) { uscoreFreq[a][it->first] /= T.size(); ucumul+= it2->second;  }
 454                                         else { uscoreFreq[a][it->first] = 0.0000; } //no user trees with that score
 455                                         //make uCumul map
 456                                         uCumul[a][it->first] = ucumul;
 457                                 }
 458
 459                                 //make rscoreFreq map and rCumul
 460                                 map<int,double>::iterator it2 = rscoreFreq[a].find(it->first);
 461                                 //get percentage of random trees with that info
 462                                 if (it2 != rscoreFreq[a].end()) {  rscoreFreq[a][it->first] /= iters; rcumul+= it2->second;  }
 463                                 else { rscoreFreq[a][it->first] = 0.0000; } //no random trees with that score
 464                                 rCumul[a][it->first] = rcumul;
 465                         }
 466
 467                         //find the signifigance of each user trees score when compared to the random trees and save for printing the summary file
 468                         for (int h = 0; h < userTreeScores[a].size(); h++) {
 469                                 UScoreSig[a].push_back(rCumul[a][userTreeScores[a][h]]);
 470                         }
 471                 }
 472
 473                 if (m->control_pressed) {
 474                                 delete reading; delete pars; delete output;
 475                                 delete tmap; for (int i = 0; i < T.size(); i++) { delete T[i]; }
 476                                 if (randomtree == "") {  outSum.close();  }
 477                                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]); } outputTypes.clear();
 478                                 m->clearGroups();
 479                                 return 0;
 480                 }
 481
 482                 //finish progress bar
 483                 reading->finish();
 484                 delete reading;
 485
 486
 487                 printParsimonyFile();
 488                 if (randomtree == "") { printUSummaryFile(); }
 489
 490                 //reset groups parameter
 491                 m->clearGroups();
 492
 493                 delete pars; delete output;
 494                 delete tmap; for (int i = 0; i < T.size(); i++) { delete T[i]; }
 495
 496                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]); } outputTypes.clear(); return 0;}
 497
 498                 m->mothurOutEndLine();
 499                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 500                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 501                 m->mothurOutEndLine();
 502
 503
 504                 return 0;
 505
 506         }
 507         catch(exception& e) {
 508                 m->errorOut(e, "ParsimonyCommand", "execute");
 509                 exit(1);
 510         }
 511 }
 512
 513 /***********************************************************/
 514 void ParsimonyCommand::printParsimonyFile() {
 515         try {
 516                 vector<double> data;
 517                 vector<string> tags;
 518
 519                 if (randomtree == "") {
 520                         tags.push_back("Score"); tags.push_back("UserFreq"); tags.push_back("UserCumul"); tags.push_back("RandFreq"); tags.push_back("RandCumul");
 521                 }else {
 522                         tags.push_back("Score"); tags.push_back("RandFreq"); tags.push_back("RandCumul");
 523                 }
 524
 525                 for(int a = 0; a < numComp; a++) {
 526                         output->initFile(groupComb[a], tags);
 527                         //print each line
 528                         for (map<int,double>::iterator it = validScores.begin(); it != validScores.end(); it++) {
 529                                 if (randomtree == "") {
 530                                         data.push_back(it->first);  data.push_back(uscoreFreq[a][it->first]); data.push_back(uCumul[a][it->first]); data.push_back(rscoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]);
 531                                 }else{
 532                                         data.push_back(it->first);  data.push_back(rscoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]);
 533                                 }
 534                                 output->output(data);
 535                                 data.clear();
 536                         }
 537                         output->resetFile();
 538                 }
 539         }
 540         catch(exception& e) {
 541                 m->errorOut(e, "ParsimonyCommand", "printParsimonyFile");
 542                 exit(1);
 543         }
 544 }
 545 /***********************************************************/
 546 int ParsimonyCommand::printUSummaryFile() {
 547         try {
 548                 //column headers
 549                 outSum << "Tree#" << '\t' << "Groups" << '\t'  <<  "ParsScore" << '\t' << "ParsSig" <<  endl;
 550                 m->mothurOut("Tree#\tGroups\tParsScore\tParsSig"); m->mothurOutEndLine();
 551
 552                 //format output
 553                 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
 554
 555
 556                 //print each line
 557                 for (int i = 0; i< T.size(); i++) {
 558                         for(int a = 0; a < numComp; a++) {
 559                                 if (m->control_pressed) {  outSum.close(); return 0; }
 560                                 if (UScoreSig[a][i] > (1/(float)iters)) {
 561                                         outSum << setprecision(6) << i+1 << '\t' << groupComb[a]  << '\t' << userTreeScores[a][i] << setprecision(itersString.length()) << '\t' << UScoreSig[a][i] << endl;
 562                                         cout << setprecision(6) << i+1 << '\t' << groupComb[a]  << '\t' << userTreeScores[a][i] << setprecision(itersString.length()) << '\t' << UScoreSig[a][i] << endl;
 563                                         m->mothurOutJustToLog(toString(i+1) + "\t" + groupComb[a] + "\t" + toString(userTreeScores[a][i]) + "\t" + toString(UScoreSig[a][i])); m->mothurOutEndLine();
 564                                 }else {
 565                                         outSum << setprecision(6) << i+1 << '\t' << groupComb[a] << '\t' << userTreeScores[a][i] << setprecision(itersString.length())  << '\t' << "<" << (1/float(iters)) << endl;
 566                                         cout << setprecision(6) << i+1 << '\t' << groupComb[a] << '\t' << userTreeScores[a][i] << setprecision(itersString.length()) << '\t' << "<" << (1/float(iters)) << endl;
 567                                         m->mothurOutJustToLog(toString(i+1) + "\t" + groupComb[a] + "\t" + toString(userTreeScores[a][i]) + "\t" + toString((1/float(iters)))); m->mothurOutEndLine();
 568                                 }
 569                         }
 570                 }
 571
 572                 outSum.close();
 573                 return 0;
 574         }
 575         catch(exception& e) {
 576                 m->errorOut(e, "ParsimonyCommand", "printUSummaryFile");
 577                 exit(1);
 578         }
 579 }
 580
 581 /***********************************************************/
 582 void ParsimonyCommand::getUserInput() {
 583         try {
 584
 585                 //create treemap
 586                 tmap = new TreeMap();
 587
 588                 m->mothurOut("Please enter the number of groups you would like to analyze: ");
 589                 cin >> numGroups;
 590                 m->mothurOutJustToLog(toString(numGroups)); m->mothurOutEndLine();
 591
 592                 int num, count;
 593                 count = 1;
 594                 numEachGroup.resize(numGroups, 0);
 595
 596
 597                 for (int i = 1; i <= numGroups; i++) {
 598                         m->mothurOut("Please enter the number of sequences in group " + toString(i) +  ": ");
 599                         cin >> num;
 600                         m->mothurOutJustToLog(toString(num)); m->mothurOutEndLine();
 601
 602                         //set tmaps seqsPerGroup
 603                         tmap->seqsPerGroup[toString(i)] = num;
 604                         tmap->addGroup(toString(i));
 605
 606                         //set tmaps namesOfSeqs
 607                         for (int j = 0; j < num; j++) {
 608                                 tmap->namesOfSeqs.push_back(toString(count));
 609                                 tmap->treemap[toString(count)].groupname = toString(i);
 610                                 count++;
 611                         }
 612                 }
 613
 614                 //clears buffer so next command doesn't have error
 615                 string s;
 616                 getline(cin, s);
 617
 618                 m->Treenames = tmap->namesOfSeqs;
 619
 620         }
 621         catch(exception& e) {
 622                 m->errorOut(e, "ParsimonyCommand", "getUserInput");
 623                 exit(1);
 624         }
 625 }
 626 /*****************************************************************/
 627 int ParsimonyCommand::readNamesFile() {
 628         try {
 629                 m->names.clear();
 630                 numUniquesInName = 0;
 631
 632                 ifstream in;
 633                 m->openInputFile(namefile, in);
 634
 635                 string first, second;
 636                 map<string, string>::iterator itNames;
 637
 638                 while(!in.eof()) {
 639                         in >> first >> second; m->gobble(in);
 640
 641                         numUniquesInName++;
 642
 643                         itNames = m->names.find(first);
 644                         if (itNames == m->names.end()) {
 645                                 m->names[first] = second;
 646
 647                                 //we need a list of names in your namefile to use above when removing extra seqs above so we don't remove them
 648                                 vector<string> dupNames;
 649                                 m->splitAtComma(second, dupNames);
 650
 651                                 for (int i = 0; i < dupNames.size(); i++) {
 652                                         nameMap[dupNames[i]] = dupNames[i];
 653                                         if ((groupfile == "") && (i != 0)) { tmap->addSeq(dupNames[i], "Group1"); }
 654                                 }
 655                         }else {  m->mothurOut(first + " has already been seen in namefile, disregarding names file."); m->mothurOutEndLine(); in.close(); m->names.clear(); namefile = ""; return 1; }
 656                 }
 657                 in.close();
 658
 659                 return 0;
 660         }
 661         catch(exception& e) {
 662                 m->errorOut(e, "ParsimonyCommand", "readNamesFile");
 663                 exit(1);
 664         }
 665 }
 666 /***********************************************************/
 667
 668