mothur.h

   1 #ifndef MOTHUR_H
   2 #define MOTHUR_H
   3
   4
   5
   6 /*
   7  *  mothur.h
   8  *  Mothur
   9  *
  10  *  Created by Sarah Westcott on 2/19/09.
  11  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
  12  *
  13  */
  14
  15 /* This file contains all the standard incudes we use in the project as well as some common utilities. */
  16
  17 //#include <cstddef>
  18
  19 //io libraries
  20 #include <iostream>
  21 #include <iomanip>
  22 #include <fstream>
  23 #include <sstream>
  24 #include <signal.h>
  25
  26 //exception
  27 #include <stdexcept>
  28 #include <exception>
  29 #include <cstdlib>
  30
  31
  32 //containers
  33 #include <vector>
  34 #include <set>
  35 #include <map>
  36 #include <string>
  37 #include <list>
  38
  39 //math
  40 #include <cmath>
  41 #include <math.h>
  42 #include <algorithm>
  43
  44 //misc
  45 #include <cerrno>
  46 #include <ctime>
  47 #include <limits>
  48 #include <readline/readline.h>
  49 #include <readline/history.h>
  50
  51 /***********************************************************************/
  52
  53 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
  54         #include <sys/wait.h>
  55         #include <unistd.h>
  56 #else
  57         #include <conio.h> //allows unbuffered screen capture from stdin
  58 #endif
  59
  60 using namespace std;
  61
  62 #define exp(x) (exp((double) x))
  63 #define sqrt(x) (sqrt((double) x))
  64 #define log10(x) (log10((double) x))
  65 #define log2(x) (log10(x)/log10(2))
  66 #define isnan(x) ((x) != (x))
  67 #define isinf(x) (fabs(x) == std::numeric_limits<double>::infinity())
  68
  69
  70 typedef unsigned long ull;
  71
  72 struct IntNode {
  73         int lvalue;
  74         int rvalue;
  75         int lcoef;
  76         int rcoef;
  77         IntNode* left;
  78         IntNode* right;
  79
  80         IntNode(int lv, int rv, IntNode* l, IntNode* r) : lvalue(lv), rvalue(rv), left(l), right(r) {};
  81         IntNode() {};
  82 };
  83
  84 struct ThreadNode {
  85         int* pid;
  86         IntNode* left;
  87         IntNode* right;
  88 };
  89
  90 /************************************************************/
  91 struct clusterNode {
  92         int numSeq;
  93         int parent;
  94         int smallChild; //used to make linkTable work with list and rabund. represents bin number of this cluster node
  95         clusterNode(int num, int par, int kid) : numSeq(num), parent(par), smallChild(kid) {};
  96 };
  97
  98 /***********************************************************************/
  99
 100 // snagged from http://www.parashift.com/c++-faq-lite/misc-technical-issues.html#faq-39.2
 101 // works for now, but there should be a way to do it without killing the whole program
 102
 103 class BadConversion : public runtime_error {
 104 public:
 105         BadConversion(const string& s) : runtime_error(s){ }
 106 };
 107
 108 //**********************************************************************************************************************
 109
 110 template<typename T>
 111 inline void convert(const string& s, T& x, bool failIfLeftoverChars = true){
 112         istringstream i(s);
 113         char c;
 114         if (!(i >> x) || (failIfLeftoverChars && i.get(c)))
 115                 throw BadConversion(s);
 116 }
 117
 118 //**********************************************************************************************************************
 119
 120 template<typename T>
 121 inline bool convertTestFloat(const string& s, T& x, bool failIfLeftoverChars = true){
 122         istringstream i(s);
 123         char c;
 124         if (!(i >> x) || (failIfLeftoverChars && i.get(c)))
 125         {
 126                 return false;
 127         }
 128         return true;
 129 }
 130
 131 //**********************************************************************************************************************
 132
 133 template<typename T>
 134 inline bool convertTest(const string& s, T& x, bool failIfLeftoverChars = true){
 135         istringstream i(s);
 136         char c;
 137         if (!(i >> x) || (failIfLeftoverChars && i.get(c)))
 138         {
 139                 cout << "unable to be converted into an integer.\n" << endl;
 140                 return false;
 141         }
 142         return true;
 143 }
 144
 145 //**********************************************************************************************************************
 146
 147 template<typename T>
 148 string toString(const T&x){
 149     stringstream output;
 150     output << x;
 151     return output.str();
 152 }
 153
 154 //**********************************************************************************************************************
 155
 156 template<typename T>
 157 string toHex(const T&x){
 158         stringstream output;
 159
 160         output << hex << x;
 161
 162     return output.str();
 163 }
 164 //**********************************************************************************************************************
 165
 166 template<typename T>
 167 string toString(const T&x, int i){
 168         stringstream output;
 169
 170         output.precision(i);
 171     output << fixed << x;
 172
 173     return output.str();
 174 }
 175 /***********************************************************************/
 176
 177 inline int openOutputFileAppend(string fileName, ofstream& fileHandle){
 178
 179         fileHandle.open(fileName.c_str(), ios::app);
 180         if(!fileHandle) {
 181                 cout << "Error: Could not open " <<  fileName << endl;
 182                 return 1;
 183         }
 184         else {
 185                 return 0;
 186         }
 187
 188 }
 189 /***********************************************************************/
 190
 191 inline void gobble(istream& f){
 192
 193         char d;
 194     while(isspace(d=f.get()))           {;}
 195         f.putback(d);
 196
 197 }
 198
 199 /***********************************************************************/
 200
 201 inline string getline(ifstream& fileHandle) {
 202         try {
 203
 204                 string line = "";
 205
 206                 while (!fileHandle.eof())       {
 207                         //get next character
 208                         char c = fileHandle.get();
 209
 210                         //are you at the end of the line
 211                         if ((c == '\n') || (c == '\r') || (c == '\f')){  break; }
 212                         else {          line += c;              }
 213                 }
 214
 215                 return line;
 216
 217         }
 218         catch(exception& e) {
 219                 cout << "Error in mothur function getline" << endl;
 220                 exit(1);
 221         }
 222 }
 223
 224 /**************************************************************************************************/
 225
 226 inline void mothurOut(string message) {
 227         try{
 228                 ofstream out;
 229                 string logFileName = "mothur.logFile";
 230                 openOutputFileAppend(logFileName, out);
 231
 232                 cout << message;
 233                 out << message;
 234
 235                 out.close();
 236         }
 237         catch(exception& e) {
 238                 cout << "Error in mothur class mothurOut" << endl;
 239                 exit(1);
 240         }
 241 }
 242 /**************************************************************************************************/
 243
 244 inline void mothurOut(string message, string precision) {
 245         try{
 246                 ofstream out;
 247                 string logFileName = "mothur.logFile";
 248                 openOutputFileAppend(logFileName, out);
 249
 250                 cout << precision << message;
 251                 out << precision << message;
 252
 253                 out.close();
 254         }
 255         catch(exception& e) {
 256                 cout << "Error in mothur class mothurOut" << endl;
 257                 exit(1);
 258         }
 259 }
 260
 261 /**************************************************************************************************/
 262
 263 inline void mothurOutEndLine() {
 264         try {
 265                 ofstream out;
 266                 string logFileName = "mothur.logFile";
 267                 openOutputFileAppend(logFileName, out);
 268
 269                 cout << endl;
 270                 out << endl;
 271
 272                 out.close();
 273         }
 274         catch(exception& e) {
 275                 cout << "error in mothur mothurOutEndLine" << endl;
 276                 exit(1);
 277         }
 278 }
 279
 280
 281 /**************************************************************************************************/
 282
 283 inline void errorOut(exception& e, string object, string function) {
 284
 285                 mothurOut("Error: ");
 286                 mothurOut(toString(e.what()));
 287                 mothurOut(" has occurred in the " + object + " class function " + function + ". Please contact Pat Schloss at mothur.bugs@gmail.com, and be sure to include the mothur.logFile with your inquiry.");
 288                 mothurOutEndLine();
 289
 290 }
 291
 292 /***********************************************************************/
 293
 294 inline bool isTrue(string f){
 295
 296         if ((f == "TRUE") || (f == "T") || (f == "true") || (f == "t")) {       return true;    }
 297         else {  return false;  }
 298 }
 299
 300 /***********************************************************************/
 301
 302 inline float roundDist(float dist, int precision){
 303
 304         return int(dist * precision + 0.5)/float(precision);
 305
 306 }
 307
 308 /***********************************************************************/
 309
 310 inline int getNumNames(string names){
 311
 312         int count = 0;
 313
 314         if(names != ""){
 315                 count = 1;
 316                 for(int i=0;i<names.size();i++){
 317                         if(names[i] == ','){
 318                                 count++;
 319                         }
 320                 }
 321         }
 322
 323         return count;
 324
 325 }
 326
 327 /**************************************************************************************************/
 328
 329 inline vector<vector<double> > binomial(int maxOrder){
 330
 331         vector<vector<double> > binomial(maxOrder+1);
 332
 333     for(int i=0;i<=maxOrder;i++){
 334                 binomial[i].resize(maxOrder+1);
 335                 binomial[i][0]=1;
 336                 binomial[0][i]=0;
 337     }
 338     binomial[0][0]=1;
 339
 340     binomial[1][0]=1;
 341     binomial[1][1]=1;
 342
 343     for(int i=2;i<=maxOrder;i++){
 344                 binomial[1][i]=0;
 345     }
 346
 347     for(int i=2;i<=maxOrder;i++){
 348                 for(int j=1;j<=maxOrder;j++){
 349                         if(i==j){       binomial[i][j]=1;                                                                       }
 350                         if(j>i) {       binomial[i][j]=0;                                                                       }
 351                         else    {       binomial[i][j]=binomial[i-1][j-1]+binomial[i-1][j];     }
 352                 }
 353     }
 354
 355         return binomial;
 356 }
 357
 358 /***********************************************************************/
 359
 360 inline string getRootName(string longName){
 361
 362         string rootName = longName;
 363
 364         if(longName.find_last_of(".") != longName.npos){
 365                 int pos = longName.find_last_of('.')+1;
 366                 rootName = longName.substr(0, pos);
 367         }
 368
 369         return rootName;
 370 }
 371 /***********************************************************************/
 372
 373 inline string getSimpleName(string longName){
 374
 375         string simpleName = longName;
 376
 377         if(longName.find_last_of("/") != longName.npos){
 378                 int pos = longName.find_last_of('/')+1;
 379                 simpleName = longName.substr(pos, longName.length());
 380         }
 381
 382         return simpleName;
 383 }
 384
 385 /***********************************************************************/
 386
 387 inline int factorial(int num){
 388         int total = 1;
 389
 390         for (int i = 1; i <= num; i++) {
 391                 total *= i;
 392         }
 393
 394         return total;
 395 }
 396 /**************************************************************************************************
 397
 398 double min(double x, double y)
 399 {
 400     if(x<y){    return x;    }
 401     else   {    return y;    }
 402 }
 403
 404 /***********************************************************************/
 405
 406 inline string getPathName(string longName){
 407
 408         string rootPathName = longName;
 409
 410         if(longName.find_last_of('/') != longName.npos){
 411                 int pos = longName.find_last_of('/')+1;
 412                 rootPathName = longName.substr(0, pos);
 413         }
 414
 415         return rootPathName;
 416 }
 417
 418 /***********************************************************************/
 419
 420 inline string getExtension(string longName){
 421
 422         string extension = longName;
 423
 424         if(longName.find_last_of('.') != longName.npos){
 425                 int pos = longName.find_last_of('.');
 426                 extension = longName.substr(pos, longName.length());
 427         }
 428
 429         return extension;
 430 }
 431 /***********************************************************************/
 432 inline bool isBlank(string fileName){
 433
 434         ifstream fileHandle;
 435         fileHandle.open(fileName.c_str());
 436         if(!fileHandle) {
 437                 mothurOut("Error: Could not open " + fileName);  mothurOutEndLine();
 438                 return false;
 439         }else {
 440                 //check for blank file
 441                 gobble(fileHandle);
 442                 if (fileHandle.eof()) { fileHandle.close(); return true;  }
 443         }
 444         return false;
 445 }
 446 /***********************************************************************/
 447
 448 inline int openInputFile(string fileName, ifstream& fileHandle){
 449
 450         fileHandle.open(fileName.c_str());
 451         if(!fileHandle) {
 452                 mothurOut("Error: Could not open " + fileName);  mothurOutEndLine();
 453                 return 1;
 454         }
 455         else {
 456                 //check for blank file
 457                 gobble(fileHandle);
 458                 if (fileHandle.eof()) { mothurOut(fileName + " is blank. Please correct."); mothurOutEndLine();  return 1;  }
 459
 460                 return 0;
 461         }
 462
 463 }
 464
 465 /***********************************************************************/
 466
 467 inline int openOutputFile(string fileName, ofstream& fileHandle){
 468
 469         fileHandle.open(fileName.c_str(), ios::trunc);
 470         if(!fileHandle) {
 471                 mothurOut("Error: Could not open " + fileName);  mothurOutEndLine();
 472                 return 1;
 473         }
 474         else {
 475                 return 0;
 476         }
 477
 478 }
 479
 480 /***********************************************************************/
 481
 482 inline int getNumSeqs(ifstream& file){
 483
 484         int numSeqs = count(istreambuf_iterator<char>(file),istreambuf_iterator<char>(), '>');
 485         file.seekg(0);
 486         return numSeqs;
 487
 488 }
 489 /***********************************************************************/
 490
 491 inline bool inVector(string member, vector<string> group){
 492
 493         for (int i = 0; i < group.size(); i++) {
 494                 if (group[i] == member) {  return true;         }
 495         }
 496
 497         return false;
 498 }
 499
 500 /***********************************************************************/
 501
 502 //This function parses the estimator options and puts them in a vector
 503 inline void splitAtDash(string& estim, vector<string>& container) {
 504         try {
 505                 string individual;
 506
 507                 while (estim.find_first_of('-') != -1) {
 508                         individual = estim.substr(0,estim.find_first_of('-'));
 509                         if ((estim.find_first_of('-')+1) <= estim.length()) { //checks to make sure you don't have dash at end of string
 510                                 estim = estim.substr(estim.find_first_of('-')+1, estim.length());
 511                                 container.push_back(individual);
 512                         }
 513                 }
 514                 //get last one
 515                 container.push_back(estim);
 516         }
 517         catch(exception& e) {
 518                 errorOut(e, "mothur", "splitAtDash");
 519                 exit(1);
 520         }
 521 }
 522
 523 /***********************************************************************/
 524 //This function parses the label options and puts them in a set
 525 inline void splitAtDash(string& estim, set<string>& container) {
 526         try {
 527                 string individual;
 528
 529                 while (estim.find_first_of('-') != -1) {
 530                         individual = estim.substr(0,estim.find_first_of('-'));
 531                         if ((estim.find_first_of('-')+1) <= estim.length()) { //checks to make sure you don't have dash at end of string
 532                                 estim = estim.substr(estim.find_first_of('-')+1, estim.length());
 533                                 container.insert(individual);
 534                         }
 535                 }
 536                 //get last one
 537                 container.insert(estim);
 538         }
 539         catch(exception& e) {
 540                 errorOut(e, "mothur", "splitAtDash");
 541                 exit(1);
 542         }
 543 }
 544 /***********************************************************************/
 545 //This function parses the line options and puts them in a set
 546 inline void splitAtDash(string& estim, set<int>& container) {
 547         try {
 548                 string individual;
 549                 int lineNum;
 550
 551                 while (estim.find_first_of('-') != -1) {
 552                         individual = estim.substr(0,estim.find_first_of('-'));
 553                         if ((estim.find_first_of('-')+1) <= estim.length()) { //checks to make sure you don't have dash at end of string
 554                                 estim = estim.substr(estim.find_first_of('-')+1, estim.length());
 555                                 convert(individual, lineNum); //convert the string to int
 556                                 container.insert(lineNum);
 557                         }
 558                 }
 559                 //get last one
 560                 convert(estim, lineNum); //convert the string to int
 561                 container.insert(lineNum);
 562         }
 563         catch(exception& e) {
 564                 errorOut(e, "mothur", "splitAtDash");
 565                 exit(1);
 566         }
 567 }
 568 /***********************************************************************/
 569 //This function parses the a string and puts peices in a vector
 570 inline void splitAtComma(string& estim, vector<string>& container) {
 571         try {
 572                 string individual;
 573
 574                 while (estim.find_first_of(',') != -1) {
 575                         individual = estim.substr(0,estim.find_first_of(','));
 576                         if ((estim.find_first_of(',')+1) <= estim.length()) { //checks to make sure you don't have comma at end of string
 577                                 estim = estim.substr(estim.find_first_of(',')+1, estim.length());
 578                                 container.push_back(individual);
 579                         }
 580                 }
 581                 //get last one
 582                 container.push_back(estim);
 583         }
 584         catch(exception& e) {
 585                 errorOut(e, "mothur", "splitAtComma");
 586                 exit(1);
 587         }
 588 }
 589 /***********************************************************************/
 590
 591 //This function splits up the various option parameters
 592 inline void splitAtComma(string& prefix, string& suffix){
 593         try {
 594                 prefix = suffix.substr(0,suffix.find_first_of(','));
 595                 if ((suffix.find_first_of(',')+2) <= suffix.length()) {  //checks to make sure you don't have comma at end of string
 596                         suffix = suffix.substr(suffix.find_first_of(',')+1, suffix.length());
 597                         string space = " ";
 598                         while(suffix.at(0) == ' ')
 599                                 suffix = suffix.substr(1, suffix.length());
 600                 }
 601
 602         }
 603         catch(exception& e) {
 604                 errorOut(e, "mothur", "splitAtComma");
 605                 exit(1);
 606         }
 607 }
 608 /***********************************************************************/
 609
 610 //This function separates the key value from the option value i.e. dist=96_...
 611 inline void splitAtEquals(string& key, string& value){
 612         try {
 613                 if(value.find_first_of('=') != -1){
 614                         key = value.substr(0,value.find_first_of('='));
 615                         if ((value.find_first_of('=')+1) <= value.length()) {
 616                                 value = value.substr(value.find_first_of('=')+1, value.length());
 617                         }
 618                 }else{
 619                         key = value;
 620                         value = 1;
 621                 }
 622         }
 623         catch(exception& e) {
 624                 errorOut(e, "mothur", "splitAtEquals");
 625                 exit(1);
 626         }
 627 }
 628 /**************************************************************************************************/
 629
 630 inline bool inUsersGroups(string groupname, vector<string> Groups) {
 631         try {
 632                 for (int i = 0; i < Groups.size(); i++) {
 633                         if (groupname == Groups[i]) { return true; }
 634                 }
 635                 return false;
 636         }
 637         catch(exception& e) {
 638                 errorOut(e, "mothur", "inUsersGroups");
 639                 exit(1);
 640         }
 641 }
 642
 643 /**************************************************************************************************/
 644
 645 inline void mothurOutJustToLog(string message) {
 646         try {
 647                 ofstream out;
 648                 string logFileName = "mothur.logFile";
 649                 openOutputFileAppend(logFileName, out);
 650
 651                 out << message;
 652
 653                 out.close();
 654         }
 655         catch(exception& e) {
 656                 errorOut(e, "mothur", "mothurOutJustToLog");
 657                 exit(1);
 658         }
 659 }
 660
 661
 662 /**************************************************************************************************/
 663
 664 inline void mothurOut(float num) {
 665         try {
 666                 ofstream out;
 667                 string logFileName = "mothur.logFile";
 668                 openOutputFileAppend(logFileName, out);
 669
 670                 cout << num;
 671                 out << num;
 672
 673                 out.close();
 674         }
 675         catch(exception& e) {
 676                 cout << "Error in mothur class mothurOut float" << endl;
 677                 exit(1);
 678         }
 679 }
 680 /***********************************************************************/
 681 inline void mothurOut(double value) {
 682         try {
 683                 ofstream out;
 684                 string logFileName = "mothur.logFile";
 685                 openOutputFileAppend(logFileName, out);
 686
 687                 cout << value;
 688                 out << value;
 689
 690                 out.close();
 691         }
 692         catch(exception& e) {
 693                 cout << "Error in mothur class mothurOut double" << endl;
 694                 exit(1);
 695         }
 696 }
 697
 698 /***********************************************************************/
 699 //this function determines if the user has given us labels that are smaller than the given label.
 700 //if so then it returns true so that the calling function can run the previous valid distance.
 701 //it's a "smart" distance function.  It also checks for invalid labels.
 702 inline bool anyLabelsToProcess(string label, set<string>& userLabels, string errorOff) {
 703         try {
 704                 set<string>::iterator it;
 705                 vector<float> orderFloat;
 706                 map<string, float> userMap;  //the conversion process removes trailing 0's which we need to put back
 707                 map<string, float>::iterator it2;
 708                 float labelFloat;
 709                 bool smaller = false;
 710
 711                 //unique is the smallest line
 712                 if (label == "unique") {  return false;  }
 713                 else { convert(label, labelFloat); }
 714
 715                 //go through users set and make them floats
 716                 for(it = userLabels.begin(); it != userLabels.end(); ++it) {
 717
 718                         float temp;
 719                         if ((*it != "unique") && (convertTestFloat(*it, temp) == true)){
 720                                 convert(*it, temp);
 721                                 orderFloat.push_back(temp);
 722                                 userMap[*it] = temp;
 723                         }else if (*it == "unique") {
 724                                 orderFloat.push_back(-1.0);
 725                                 userMap["unique"] = -1.0;
 726                         }else {
 727                                 if (errorOff == "") {  mothurOut(*it + " is not a valid label."); mothurOutEndLine();  }
 728                                 userLabels.erase(*it);
 729                                 it--;
 730                         }
 731                 }
 732
 733                 //sort order
 734                 sort(orderFloat.begin(), orderFloat.end());
 735
 736                 /*************************************************/
 737                 //is this label bigger than any of the users labels
 738                 /*************************************************/
 739
 740                 //loop through order until you find a label greater than label
 741                 for (int i = 0; i < orderFloat.size(); i++) {
 742                         if (orderFloat[i] < labelFloat) {
 743                                 smaller = true;
 744                                 if (orderFloat[i] == -1) {
 745                                         if (errorOff == "") { mothurOut("Your file does not include the label unique."); mothurOutEndLine(); }
 746                                         userLabels.erase("unique");
 747                                 }
 748                                 else {
 749                                         if (errorOff == "") { mothurOut("Your file does not include the label "); mothurOutEndLine(); }
 750                                         string s = "";
 751                                         for (it2 = userMap.begin(); it2!= userMap.end(); it2++) {
 752                                                 if (it2->second == orderFloat[i]) {
 753                                                         s = it2->first;
 754                                                         //remove small labels
 755                                                         userLabels.erase(s);
 756                                                         break;
 757                                                 }
 758                                         }
 759                                         if (errorOff == "") { mothurOut(s + ". I will use the next smallest distance. "); mothurOutEndLine(); }
 760                                 }
 761                         //since they are sorted once you find a bigger one stop looking
 762                         }else { break; }
 763                 }
 764
 765                 return smaller;
 766
 767         }
 768         catch(exception& e) {
 769                 errorOut(e, "mothur", "anyLabelsToProcess");
 770                 exit(1);
 771         }
 772 }
 773
 774 /**************************************************************************************************/
 775 #endif
 776