]> git.donarmstrong.com Git - mothur.git/blob - distancecommand.cpp
added MPI code, broke up chimera.seqs into 5 separated commands, added parse.sff...
[mothur.git] / distancecommand.cpp
1 /*
2  *  distancecommand.cpp
3  *  Mothur
4  *
5  *  Created by Sarah Westcott on 5/7/09.
6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
7  *
8  */
9
10 #include "distancecommand.h"
11 #include "ignoregaps.h"
12 #include "eachgapdist.h"
13 #include "eachgapignore.h"
14 #include "onegapdist.h"
15 #include "onegapignore.h"
16
17 //**********************************************************************************************************************
18
19 DistanceCommand::DistanceCommand(string option) {
20         try {
21                 abort = false;
22                 Estimators.clear();
23                                 
24                 //allow user to run help
25                 if(option == "help") { help(); abort = true; }
26                 
27                 else {
28                         //valid paramters for this command
29                         string Array[] =  {"fasta", "output", "calc", "countends", "cutoff", "processors", "outputdir","inputdir"};
30                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
31                         
32                         OptionParser parser(option);
33                         map<string, string> parameters = parser.getParameters();
34                         
35                         ValidParameters validParameter;
36                         map<string, string>::iterator it2;
37                 
38                         //check to make sure all parameters are valid for command
39                         for (it2 = parameters.begin(); it2 != parameters.end(); it2++) { 
40                                 if (validParameter.isValidParameter(it2->first, myArray, it2->second) != true) {  abort = true;  }
41                         }
42                         
43                         //if the user changes the input directory command factory will send this info to us in the output parameter 
44                         string inputDir = validParameter.validFile(parameters, "inputdir", false);              
45                         if (inputDir == "not found"){   inputDir = "";          }
46                         else {
47                                 string path;
48                                 it2 = parameters.find("fasta");
49                                 //user has given a template file
50                                 if(it2 != parameters.end()){ 
51                                         path = hasPath(it2->second);
52                                         //if the user has not given a path then, add inputdir. else leave path alone.
53                                         if (path == "") {       parameters["fasta"] = inputDir + it2->second;           }
54                                 }
55                         }
56
57                         //check for required parameters
58                         fastafile = validParameter.validFile(parameters, "fasta", true);
59                         if (fastafile == "not found") { m->mothurOut("fasta is a required parameter for the dist.seqs command."); m->mothurOutEndLine(); abort = true; }
60                         else if (fastafile == "not open") { abort = true; }     
61                         else{
62                                 ifstream inFASTA;
63                                 openInputFile(fastafile, inFASTA);
64                                 alignDB = SequenceDB(inFASTA); 
65                                 inFASTA.close();
66                         }
67                         
68                         //if the user changes the output directory command factory will send this info to us in the output parameter 
69                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  
70                                 outputDir = ""; 
71                                 outputDir += hasPath(fastafile); //if user entered a file with a path then preserve it  
72                         }
73
74                         //check for optional parameter and set defaults
75                         // ...at some point should added some additional type checking...
76                         calc = validParameter.validFile(parameters, "calc", false);                     
77                         if (calc == "not found") { calc = "onegap";  }
78                         else { 
79                                  if (calc == "default")  {  calc = "onegap";  }
80                         }
81                         splitAtDash(calc, Estimators);
82
83                         string temp;
84                         temp = validParameter.validFile(parameters, "countends", false);        if(temp == "not found"){        temp = "T";     }
85                         convert(temp, countends); 
86                         
87                         temp = validParameter.validFile(parameters, "cutoff", false);           if(temp == "not found"){        temp = "1.0"; }
88                         convert(temp, cutoff); 
89                         
90                         temp = validParameter.validFile(parameters, "processors", false);       if(temp == "not found"){        temp = "1"; }
91                         convert(temp, processors); 
92                         
93                         output = validParameter.validFile(parameters, "output", false);         if(output == "not found"){      output = "column"; }
94                         
95                         if ((output != "column") && (output != "lt") && (output != "square")) { m->mothurOut(output + " is not a valid output form. Options are column, lt and square. I will use column."); m->mothurOutEndLine(); output = "column"; }
96                         
97                         ValidCalculators validCalculator;
98                         
99                         if (isTrue(countends) == true) {
100                                 for (int i=0; i<Estimators.size(); i++) {
101                                         if (validCalculator.isValidCalculator("distance", Estimators[i]) == true) { 
102                                                 if (Estimators[i] == "nogaps")                  {       distCalculator = new ignoreGaps();      }
103                                                 else if (Estimators[i] == "eachgap")    {       distCalculator = new eachGapDist();     }
104                                                 else if (Estimators[i] == "onegap")             {       distCalculator = new oneGapDist();      }
105                                         }
106                                 }
107                         }else {
108                                 for (int i=0; i<Estimators.size(); i++) {
109                                         if (validCalculator.isValidCalculator("distance", Estimators[i]) == true) { 
110                                                 if (Estimators[i] == "nogaps")          {       distCalculator = new ignoreGaps();                                      }
111                                                 else if (Estimators[i] == "eachgap"){   distCalculator = new eachGapIgnoreTermGapDist();        }
112                                                 else if (Estimators[i] == "onegap")     {       distCalculator = new oneGapIgnoreTermGapDist();         }
113                                         }
114                                 }
115                         }
116
117                 }
118                                 
119         }
120         catch(exception& e) {
121                 m->errorOut(e, "DistanceCommand", "DistanceCommand");
122                 exit(1);
123         }
124 }
125
126 //**********************************************************************************************************************
127
128 DistanceCommand::~DistanceCommand(){
129         
130         for(int i=0;i<lines.size();i++){
131                 delete lines[i];
132         }
133         
134 }
135         
136 //**********************************************************************************************************************
137
138 void DistanceCommand::help(){
139         try {
140                 m->mothurOut("The dist.seqs command reads a file containing sequences and creates a distance file.\n");
141                 m->mothurOut("The dist.seqs command parameters are fasta, calc, countends, output, cutoff and processors.  \n");
142                 m->mothurOut("The fasta parameter is required.\n");
143                 m->mothurOut("The calc parameter allows you to specify the method of calculating the distances.  Your options are: nogaps, onegap or eachgap. The default is onegap.\n");
144                 m->mothurOut("The countends parameter allows you to specify whether to include terminal gaps in distance.  Your options are: T or F. The default is T.\n");
145                 m->mothurOut("The cutoff parameter allows you to specify maximum distance to keep. The default is 1.0.\n");
146                 m->mothurOut("The output parameter allows you to specify format of your distance matrix. Options are column, lt, and square. The default is column.\n");
147                 m->mothurOut("The processors parameter allows you to specify number of processors to use.  The default is 1.\n");
148                 m->mothurOut("The dist.seqs command should be in the following format: \n");
149                 m->mothurOut("dist.seqs(fasta=yourFastaFile, calc=yourCalc, countends=yourEnds, cutoff= yourCutOff, processors=yourProcessors) \n");
150                 m->mothurOut("Example dist.seqs(fasta=amazon.fasta, calc=eachgap, countends=F, cutoff= 2.0, processors=3).\n");
151                 m->mothurOut("Note: No spaces between parameter labels (i.e. calc), '=' and parameters (i.e.yourCalc).\n\n");
152         }
153         catch(exception& e) {
154                 m->errorOut(e, "DistanceCommand", "help");
155                 exit(1);
156         }
157 }
158 //**********************************************************************************************************************
159
160 int DistanceCommand::execute(){
161         try {
162                 
163                 if (abort == true) { return 0; }
164                 
165                 int startTime = time(NULL);
166                 
167                 int numSeqs = alignDB.getNumSeqs();
168                 cutoff += 0.005;
169                 
170                 string outputFile;
171                                 
172                 if (output == "lt") { //does the user want lower triangle phylip formatted file 
173                         outputFile = outputDir + getRootName(getSimpleName(fastafile)) + "phylip.dist";
174                         remove(outputFile.c_str());
175                         
176                         //output numSeqs to phylip formatted dist file
177                 }else if (output == "column") { //user wants column format
178                         outputFile = outputDir + getRootName(getSimpleName(fastafile)) + "dist";
179                         remove(outputFile.c_str());
180                 }else { //assume square
181                         outputFile = outputDir + getRootName(getSimpleName(fastafile)) + "square.dist";
182                         remove(outputFile.c_str());
183                 }
184                 
185
186 #ifdef USE_MPI
187                 
188                 int pid, start, end; 
189                 int tag = 2001;
190                                 
191                 MPI_Status status; 
192                 MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
193                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
194                 
195                 //each process gets where it should start and stop in the file
196                 start = int (sqrt(float(pid)/float(processors)) * numSeqs);
197                 end = int (sqrt(float(pid+1)/float(processors)) * numSeqs);
198                 
199                 if (output != "lt") {
200                         MPI_File outMPI;
201                         int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY; 
202                         
203                         char filename[outputFile.length()];
204                         strcpy(filename, outputFile.c_str());
205                         
206                         MPI_File_open(MPI_COMM_WORLD, filename, amode, MPI_INFO_NULL, &outMPI);
207                         
208                         if (m->control_pressed) {   MPI_File_close(&outMPI);  delete distCalculator;  return 0;  }
209
210                         if (pid == 0) { //you are the root process 
211                         
212                                 //do your part
213                                 string outputMyPart;
214                                 driverMPI(start, end, outMPI, cutoff);
215                                 
216                                 if (m->control_pressed) { MPI_File_close(&outMPI);  delete distCalculator;  return 0; }
217                         
218                                 //wait on chidren
219                                 for(int i = 1; i < processors; i++) { 
220                                         if (m->control_pressed) { MPI_File_close(&outMPI);  delete distCalculator;  return 0; }
221                                         
222                                         char buf[4];
223                                         MPI_Recv(buf, 4, MPI_CHAR, i, tag, MPI_COMM_WORLD, &status); 
224                                 }
225                         }else { //you are a child process
226                                 //do your part
227                                 driverMPI(start, end, outMPI, cutoff);
228                                 
229                                 if (m->control_pressed) { MPI_File_close(&outMPI);  delete distCalculator;  return 0; }
230                         
231                                 char buf[4];
232                                 strcpy(buf, "done"); 
233                                 //tell parent you are done.
234                                 MPI_Send(buf, 4, MPI_CHAR, 0, tag, MPI_COMM_WORLD);
235                         }
236                         
237                         MPI_File_close(&outMPI);
238                         
239                 }else { //lower triangle format
240                         if (pid == 0) { //you are the root process 
241                         
242                                 //do your part
243                                 string outputMyPart;
244                                 long mySize;
245                                 driverMPI(start, end, outputFile, mySize);
246         
247                                 if (m->control_pressed) {  delete distCalculator;  return 0; }
248                                 
249                                 int amode=MPI_MODE_APPEND|MPI_MODE_WRONLY|MPI_MODE_CREATE; //
250                                 MPI_File outMPI;
251                                 MPI_File inMPI;
252                         
253                                 char filename[outputFile.length()];
254                                 strcpy(filename, outputFile.c_str());
255                         
256                                 MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
257
258                                 //wait on chidren
259                                 for(int b = 1; b < processors; b++) { 
260                                         long fileSize;
261                                         
262                                         if (m->control_pressed) { MPI_File_close(&outMPI);  delete distCalculator;  return 0; }
263                                         
264                                         MPI_Recv(&fileSize, 1, MPI_LONG, b, tag, MPI_COMM_WORLD, &status); 
265                                         
266                                         string outTemp = outputFile + toString(b) + ".temp";
267                                         char buf[outTemp.length()];
268                                         strcpy(buf, outTemp.c_str());
269                                         
270                                         MPI_File_open(MPI_COMM_SELF, buf, MPI_MODE_DELETE_ON_CLOSE|MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);
271                                         
272                                         int count = 0;
273                                         while (count < fileSize) { //read 1000 characters at a time
274                                                 //send freqs
275                                                 char buf2[1];
276                                                 MPI_File_read(inMPI, buf2, 1, MPI_CHAR, &status);
277                                                 MPI_File_write(outMPI, buf2, 1, MPI_CHAR, &status);
278                                                 count += 1;
279                                         }
280                                         
281                                         MPI_File_close(&inMPI); //deleted on close
282                                 }
283                                 
284                                 MPI_File_close(&outMPI);
285                         }else { //you are a child process
286                                 //do your part
287                                 long size;
288                                 driverMPI(start, end, (outputFile + toString(pid) + ".temp"), size);
289                                 
290                                 if (m->control_pressed) { delete distCalculator;  return 0; }
291                         
292                                 //tell parent you are done.
293                                 MPI_Send(&size, 1, MPI_LONG, 0, tag, MPI_COMM_WORLD);
294                         }
295                 }
296 #else           
297                                 
298         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
299                 //if you don't need to fork anything
300                 if(processors == 1){
301                         driver(0, numSeqs, outputFile, cutoff);
302                 }else{ //you have multiple processors
303                         
304                         for (int i = 0; i < processors; i++) {
305                                 lines.push_back(new linePair());
306                                 lines[i]->start = int (sqrt(float(i)/float(processors)) * numSeqs);
307                                 lines[i]->end = int (sqrt(float(i+1)/float(processors)) * numSeqs);
308                         }
309
310                         createProcesses(outputFile); 
311                 
312                         map<int, int>::iterator it = processIDS.begin();
313                         rename((outputFile + toString(it->second) + ".temp").c_str(), outputFile.c_str());
314                         it++;
315                         
316                         //append and remove temp files
317                         for (; it != processIDS.end(); it++) {
318                                 appendFiles((outputFile + toString(it->second) + ".temp"), outputFile);
319                                 remove((outputFile + toString(it->second) + ".temp").c_str());
320                         }
321                 }
322         #else
323                 ifstream inFASTA;
324                 driver(0, numSeqs, outputFile, cutoff);
325         #endif
326         
327 #endif
328                 if (m->control_pressed) { delete distCalculator; remove(outputFile.c_str()); return 0; }
329                 
330                 #ifdef USE_MPI
331                         MPI_Comm_rank(MPI_COMM_WORLD, &pid); 
332                                         
333                         if (pid == 0) { //only one process should output to screen
334                 #endif
335                 
336                 if (output == "square") {  convertMatrix(outputFile); }
337                 
338                 #ifdef USE_MPI
339                         }
340                 #endif
341                 
342                 if (m->control_pressed) { delete distCalculator; remove(outputFile.c_str()); return 0; }
343                 
344                 delete distCalculator;
345                 
346                 m->mothurOutEndLine();
347                 m->mothurOut("Output File Name: "); m->mothurOutEndLine();
348                 m->mothurOut(outputFile); m->mothurOutEndLine();
349                 m->mothurOutEndLine();
350                 m->mothurOut("It took " + toString(time(NULL) - startTime) + " to calculate the distances for " + toString(numSeqs) + " sequences."); m->mothurOutEndLine();
351                 return 0;
352                 
353         }
354         catch(exception& e) {
355                 m->errorOut(e, "DistanceCommand", "execute");
356                 exit(1);
357         }
358 }
359 /**************************************************************************************************/
360 void DistanceCommand::createProcesses(string filename) {
361         try {
362 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
363                 int process = 0;
364                 processIDS.clear();
365                 
366                 //loop through and create all the processes you want
367                 while (process != processors) {
368                         int pid = fork();
369                         
370                         if (pid > 0) {
371                                 processIDS[lines[process]->end] = pid;  //create map from line number to pid so you can append files in correct order later
372                                 process++;
373                         }else if (pid == 0){
374                                 driver(lines[process]->start, lines[process]->end, filename + toString(getpid()) + ".temp", cutoff);
375                                 exit(0);
376                         }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
377                 }
378         
379                 //force parent to wait until all the processes are done
380                 for (map<int, int>::iterator it = processIDS.begin(); it != processIDS.end(); it++) { 
381                         int temp = it->second;
382                         wait(&temp);
383                 }
384 #endif
385         }
386         catch(exception& e) {
387                 m->errorOut(e, "DistanceCommand", "createProcesses");
388                 exit(1);
389         }
390 }
391
392 /**************************************************************************************************/
393 /////// need to fix to work with calcs and sequencedb
394 int DistanceCommand::driver(int startLine, int endLine, string dFileName, float cutoff){
395         try {
396
397                 int startTime = time(NULL);
398                 
399                 //column file
400                 ofstream outFile(dFileName.c_str(), ios::trunc);
401                 outFile.setf(ios::fixed, ios::showpoint);
402                 outFile << setprecision(4);
403                 
404                 if((output == "lt") && startLine == 0){ outFile << alignDB.getNumSeqs() << endl;        }
405                 
406                 for(int i=startLine;i<endLine;i++){
407                         if(output == "lt")      {       
408                                 string name = alignDB.get(i).getName();
409                                 if (name.length() < 10) { //pad with spaces to make compatible
410                                         while (name.length() < 10) {  name += " ";  }
411                                 }
412                                 outFile << name << '\t';        
413                         }
414                         for(int j=0;j<i;j++){
415                                 
416                                 if (m->control_pressed) { outFile.close(); return 0;  }
417                                 
418                                 distCalculator->calcDist(alignDB.get(i), alignDB.get(j));
419                                 double dist = distCalculator->getDist();
420                                 
421                                 if(dist <= cutoff){
422                                         if (output == "column") { outFile << alignDB.get(i).getName() << ' ' << alignDB.get(j).getName() << ' ' << dist << endl; }
423                                 }
424                                 if (output == "lt") {  outFile << dist << '\t'; }
425                                 
426                                 if (output == "square") { //make a square column you can convert to square phylip
427                                         outFile << alignDB.get(i).getName() << '\t' << alignDB.get(j).getName() << '\t' << dist << endl;
428                                         outFile << alignDB.get(j).getName() << '\t' << alignDB.get(i).getName() << '\t' << dist << endl;
429                                 }
430
431                         }
432                         
433                         if (output == "lt") { outFile << endl; }
434                         
435                         if(i % 100 == 0){
436                                 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
437                         }
438                         
439                 }
440                 m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
441                 
442                 outFile.close();
443                 
444                 return 1;
445         }
446         catch(exception& e) {
447                 m->errorOut(e, "DistanceCommand", "driver");
448                 exit(1);
449         }
450 }
451 #ifdef USE_MPI
452 /**************************************************************************************************/
453 /////// need to fix to work with calcs and sequencedb
454 int DistanceCommand::driverMPI(int startLine, int endLine, MPI_File& outMPI, float cutoff){
455         try {
456                 MPI_Status status;
457                 int startTime = time(NULL);
458                 
459                 string outputString = "";
460                 
461                 for(int i=startLine;i<endLine;i++){
462         
463                         for(int j=0;j<i;j++){
464                                 
465                                 if (m->control_pressed) {  return 0;  }
466                                 
467                                 distCalculator->calcDist(alignDB.get(i), alignDB.get(j));
468                                 double dist = distCalculator->getDist();
469                                 
470                                 if(dist <= cutoff){
471                                         if (output == "column") { outputString += (alignDB.get(i).getName() + ' ' + alignDB.get(j).getName() + ' ' + toString(dist) + '\n'); }
472                                 }
473                                 
474                                 if (output == "square") { //make a square column you can convert to square phylip
475                                         outputString += (alignDB.get(i).getName() + ' ' + alignDB.get(j).getName() + ' ' + toString(dist) + '\n');
476                                         outputString += (alignDB.get(j).getName() + ' ' + alignDB.get(i).getName() + ' ' + toString(dist) + '\n');
477                                 }
478                         }
479                         
480                         if(i % 100 == 0){
481                                 //m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
482                                 cout << i << '\t' << (time(NULL) - startTime) << endl;
483                         }
484                         
485                          
486                         //send results to parent
487                         int length = outputString.length();
488                         char buf[length];
489                         strcpy(buf, outputString.c_str()); 
490                         
491                         MPI_File_write_shared(outMPI, buf, length, MPI_CHAR, &status);
492                         outputString = "";
493                         
494                 }
495                 
496                 //m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
497                 cout << (endLine-1) << '\t' << (time(NULL) - startTime) << endl;                
498                 return 1;
499         }
500         catch(exception& e) {
501                 m->errorOut(e, "DistanceCommand", "driverMPI");
502                 exit(1);
503         }
504 }
505 /**************************************************************************************************/
506 /////// need to fix to work with calcs and sequencedb
507 int DistanceCommand::driverMPI(int startLine, int endLine, string file, long& size){
508         try {
509                 MPI_Status status;
510                 
511                 MPI_File outMPI;
512                 int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY; 
513                 
514                 char filename[file.length()];
515                 strcpy(filename, file.c_str());
516                 
517                 MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
518
519                 int startTime = time(NULL);
520                 
521                 string outputString = "";
522                 size = 0;
523                 
524                 if((output == "lt") && startLine == 0){ outputString += toString(alignDB.getNumSeqs()) + "\n";  }
525                 
526                 for(int i=startLine;i<endLine;i++){
527                         if(output == "lt")      {       
528                                 string name = alignDB.get(i).getName();
529                                 if (name.length() < 10) { //pad with spaces to make compatible
530                                         while (name.length() < 10) {  name += " ";  }
531                                 }
532                                 outputString += name + "\t";    
533                         }
534                         for(int j=0;j<i;j++){
535                                 
536                                 if (m->control_pressed) {  return 0;  }
537                                 
538                                 distCalculator->calcDist(alignDB.get(i), alignDB.get(j));
539                                 double dist = distCalculator->getDist();
540                                 
541                                 if (output == "lt") {  outputString += toString(dist) + "\t"; }
542                         }
543                         
544                         if (output == "lt") { outputString += "\n"; }
545
546                 
547                         if(i % 100 == 0){
548                                 //m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
549                                 cout << i << '\t' << (time(NULL) - startTime) << endl;
550                         }
551                         
552                         
553                         //send results to parent
554                         int length = outputString.length();
555                         char buf[length];
556                         strcpy(buf, outputString.c_str()); 
557                         
558                         MPI_File_write(outMPI, buf, length, MPI_CHAR, &status);
559                         size += outputString.length();
560                         outputString = "";
561                         
562                         
563                 }
564                 
565                 //m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
566                 cout << (endLine-1) << '\t' << (time(NULL) - startTime) << endl;
567                 MPI_File_close(&outMPI);
568                 
569                 return 1;
570         }
571         catch(exception& e) {
572                 m->errorOut(e, "DistanceCommand", "driverMPI");
573                 exit(1);
574         }
575 }
576 #endif
577 /**************************************************************************************************/
578 int DistanceCommand::convertMatrix(string outputFile) {
579         try{
580
581                 //sort file by first column so the distances for each row are together
582                 string outfile = getRootName(outputFile) + "sorted.dist.temp";
583                 
584                 //use the unix sort 
585                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
586                         string command = "sort -n " + outputFile + " -o " + outfile;
587                         system(command.c_str());
588                 #else //sort using windows sort
589                         string command = "sort " + outputFile + " /O " + outfile;
590                         system(command.c_str());
591                 #endif
592                 
593
594                 //output to new file distance for each row and save positions in file where new row begins
595                 ifstream in;
596                 openInputFile(outfile, in);
597                 
598                 ofstream out;
599                 openOutputFile(outputFile, out);
600                 
601                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
602
603                 out << alignDB.getNumSeqs() << endl;
604                 
605                 //get first currentRow
606                 string first, currentRow, second;
607                 float dist;
608                 map<string, float> rowDists; //take advantage of the fact that maps are already sorted by key 
609                 map<string, float>::iterator it;
610                 
611                 in >> first;
612                 currentRow = first;
613                 
614                 rowDists[first] = 0.00; //distance to yourself is 0.0
615                 
616                 in.seekg(0);
617                 //openInputFile(outfile, in);
618                 
619                 while(!in.eof()) {
620                         if (m->control_pressed) { in.close(); remove(outfile.c_str()); out.close(); return 0; }
621                         
622                         in >> first >> second >> dist; gobble(in);
623                                 
624                         if (first != currentRow) {
625                                 //print out last row
626                                 out << currentRow << '\t'; //print name
627
628                                 //print dists
629                                 for (it = rowDists.begin(); it != rowDists.end(); it++) {
630                                         out << it->second << '\t';
631                                 }
632                                 out << endl;
633                                 
634                                 //start new row
635                                 currentRow = first;
636                                 rowDists.clear();
637                                 rowDists[first] = 0.00;
638                                 rowDists[second] = dist;
639                         }else{
640                                 rowDists[second] = dist;
641                         }
642                 }
643                 //print out last row
644                 out << currentRow << '\t'; //print name
645                                 
646                 //print dists
647                 for (it = rowDists.begin(); it != rowDists.end(); it++) {
648                         out << it->second << '\t';
649                 }
650                 out << endl;
651                 
652                 in.close();
653                 out.close();
654                 
655                 remove(outfile.c_str());
656                 
657                 return 1;
658                 
659         }
660         catch(exception& e) {
661                 m->errorOut(e, "DistanceCommand", "convertMatrix");
662                 exit(1);
663         }
664 }
665 /**************************************************************************************************/
666 int DistanceCommand::convertToLowerTriangle(string outputFile) {
667         try{
668
669                 //sort file by first column so the distances for each row are together
670                 string outfile = getRootName(outputFile) + "sorted.dist.temp";
671                 
672                 //use the unix sort 
673                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
674                         string command = "sort -n " + outputFile + " -o " + outfile;
675                         system(command.c_str());
676                 #else //sort using windows sort
677                         string command = "sort " + outputFile + " /O " + outfile;
678                         system(command.c_str());
679                 #endif
680                 
681
682                 //output to new file distance for each row and save positions in file where new row begins
683                 ifstream in;
684                 openInputFile(outfile, in);
685                 
686                 ofstream out;
687                 openOutputFile(outputFile, out);
688                 
689                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
690
691                 out << alignDB.getNumSeqs() << endl;
692                 
693                 //get first currentRow
694                 string first, currentRow, second;
695                 float dist;
696                 int i, j;
697                 i = 0; j = 0;
698                 map<string, float> rowDists; //take advantage of the fact that maps are already sorted by key 
699                 map<string, float>::iterator it;
700                 
701                 in >> first;
702                 currentRow = first;
703                 
704                 rowDists[first] = 0.00; //distance to yourself is 0.0
705                 
706                 in.seekg(0);
707                 //openInputFile(outfile, in);
708                 
709                 while(!in.eof()) {
710                         if (m->control_pressed) { in.close(); remove(outfile.c_str()); out.close(); return 0; }
711                         
712                         in >> first >> second >> dist; gobble(in);
713                                 
714                         if (first != currentRow) {
715                                 //print out last row
716                                 out << currentRow << '\t'; //print name
717
718                                 //print dists
719                                 for (it = rowDists.begin(); it != rowDists.end(); it++) {
720                                         if (j >= i) { break; }
721                                         out << it->second << '\t';
722                                         j++;
723                                 }
724                                 out << endl;
725                                 
726                                 //start new row
727                                 currentRow = first;
728                                 rowDists.clear();
729                                 rowDists[first] = 0.00;
730                                 rowDists[second] = dist;
731                                 j = 0;
732                                 i++;
733                         }else{
734                                 rowDists[second] = dist;
735                         }
736                 }
737                 //print out last row
738                 out << currentRow << '\t'; //print name
739                                 
740                 //print dists
741                 for (it = rowDists.begin(); it != rowDists.end(); it++) {
742                         out << it->second << '\t';
743                 }
744                 out << endl;
745                 
746                 in.close();
747                 out.close();
748                 
749                 remove(outfile.c_str());
750                 
751                 return 1;
752                 
753         }
754         catch(exception& e) {
755                 m->errorOut(e, "DistanceCommand", "convertToLowerTriangle");
756                 exit(1);
757         }
758 }
759 /**************************************************************************************************
760 void DistanceCommand::appendFiles(string temp, string filename) {
761         try{
762                 ofstream output;
763                 ifstream input;
764         
765                 //open output file in append mode
766                 openOutputFileAppend(filename, output);
767                 openInputFile(temp, input);
768                 
769                 while(char c = input.get()){
770                         if(input.eof())         {       break;                  }
771                         else                            {       output << c;    }
772                 }
773                 
774                 input.close();
775                 output.close();
776         }
777         catch(exception& e) {
778                 m->errorOut(e, "DistanceCommand", "appendFiles");
779                 exit(1);
780         }
781 }
782 /**************************************************************************************************/
783
784