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