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