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