tools.h

/****************************************************************************
 *   Copyright (C) 2022 by
 *   Nhan Ly-Trong <trongnhan.uit@gmail.com>
 *   Chris Bielow <chris.bielow@fu-berlin.de>
 *   Nicola De Maio <demaio@ebi.ac.uk>
 *   BUI Quang Minh <m.bui@anu.edu.au>
 *
 *
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the
 *   Free Software Foundation, Inc.,
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 ***************************************************************************/
 
#pragma once
 
#include <cmaple_config.h>
#include <sys/stat.h>
#include <algorithm>
#include <cfloat>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <map>
#include <random>
#include <regex>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <thread>
#include <vector>
#include <cassert>
#include "operatingsystem.h"
#ifdef _OPENMP
#include <omp.h> /* OpenMP */
#endif
 
#ifdef NDEBUG
#define ASSERT(EXPRESSION) ((void)0)
#else
#if defined(__GNUC__) || defined(__clang__)
#define ASSERT(EXPRESSION)                                             \
  ((EXPRESSION) ? (void)0                                              \
                : cmaple::_my_assert(#EXPRESSION, __PRETTY_FUNCTION__, \
                                     __FILE__, __LINE__))
#else
#define ASSERT(EXPRESSION) \
  ((EXPRESSION)            \
       ? (void)0           \
       : cmaple::_my_assert(#EXPRESSION, __func__, __FILE__, __LINE__))
#endif
#endif
 
#define USE_HASH_MAP
 
#if defined(__GNUC__) && !defined(GCC_VERSION)
#define GCC_VERSION \
  (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__)
// #else
// #define GCC_VERSION 0
#endif
 
// for MSVC
#ifndef __func__
#define __func__ __FUNCTION__
#endif
 
#if defined(USE_HASH_MAP)
//    #include <unordered_map>
//    #include <unordered_set>
 
#if defined(_MSC_VER)
#include <unordered_map>
#include <unordered_set>
#elif defined(__clang__)
// libc++ detected:     _LIBCPP_VERSION
// libstdc++ detected:  __GLIBCXX__
#if __has_include(<unordered_map>)  // defines _LIBCPP_VERSION
#include <unordered_map>
#include <unordered_set>
#else
#include <tr1/unordered_map>
#include <tr1/unordered_set>
using namespace std::tr1;
#endif
#elif !defined(__GNUC__)
#include <hash_map>
#include <hash_set>
using namespace stdext;
#elif GCC_VERSION < 40300
#include <ext/hash_map>
#include <ext/hash_set>
using namespace __gnu_cxx;
#define unordered_map hash_map
#define unordered_set hash_set
/*
NHANLT: resolve ambiguous reference
std::tr1 is Technical Report 1, a proposal of extensions for C++0X when waiting
for C++11 to get approved. Most of features (including
std::tr1::unordered_map/set) in Technical Report 1 had been merged into C++11.
C++11 had been fully implemented since GCC 4.8.1 -> we only need to use std::tr1
in GCC older than 4.8.1
*/
#elif GCC_VERSION < 40800
#include <tr1/unordered_map>
#include <tr1/unordered_set>
using namespace std::tr1;
#else
#include <unordered_map>
#include <unordered_set>
#endif
 
#else
#include <map>
#include <set>
#endif
 
#if defined(USE_HASH_MAP) && GCC_VERSION < 40300 && !defined(_MSC_VER) && \
    !defined(__clang__)
/*
        Define the hash function of Split
 */
#if !defined(__GNUC__)
namespace stdext
#else
namespace __gnu_cxx
#endif
{
template <>
struct hash<string> {
  size_t operator()(string str) const {
    hash<const char*> hash_str;
    return hash_str(str.c_str());
  }
};
}  // namespace stdext
#endif  // USE_HASH_MAP
 
namespace cmaple {
// redefine assertion
inline void _my_assert(const char* expression,
                       const char* func,
                       const char* file,
                       int line) {
  char* sfile = const_cast<char*>(strrchr(file, '/'));
  if (!sfile)
    sfile = const_cast<char*>(file);
  else
    sfile++;
  std::cerr << sfile << ":" << line << ": " << func << ": Assertion `"
            << expression << "' failed." << std::endl;
  abort();
}
 
typedef uint16_t StateType;
 
typedef int32_t PositionType;
 
typedef uint32_t NumSeqsType;
 
typedef int16_t LengthType;
typedef int32_t LengthTypeLarge;
 
typedef double RealNumType;
 
typedef std::vector<RealNumType> RealNumberVector;
 
typedef std::vector<int> IntList;
 
typedef std::vector<int> IntVector;
 
typedef std::vector<bool> BoolVector;
 
typedef std::vector<char> CharVector;
 
typedef std::vector<std::string> StrVector;
 
typedef unsigned int UINT;
 
enum RegionType : StateType {
  TYPE_R = 250,
  TYPE_O,
  TYPE_N,
  TYPE_DEL,
  TYPE_INVALID
};
 
enum VerboseMode { VB_QUIET, VB_MIN, VB_MED, VB_MAX, VB_DEBUG };
 
extern VerboseMode verbose_mode;
 
enum MiniIndex : NumSeqsType {
  TOP,
  LEFT,
  RIGHT,
  UNDEFINED  // to handle case return node = null
};
 
/*--------------------------- NODE's INDEX -----------------------------------*/
struct Index {
    Index() : vector_index_(0), mini_index_(UNDEFINED){};
 
  Index(NumSeqsType vec_index, MiniIndex mini_index)
      : vector_index_(vec_index), mini_index_(mini_index){};
 
  MiniIndex getMiniIndex() const { return mini_index_; }
 
  MiniIndex getFlipMiniIndex() const { return MiniIndex(3 - mini_index_); }
 
  void setMiniIndex(MiniIndex mini_index) { mini_index_ = mini_index; }
 
  NumSeqsType getVectorIndex() const { return vector_index_; }
 
  void setVectorIndex(NumSeqsType vec_index) { vector_index_ = vec_index; }
 
  bool operator==(const Index& rhs) const {
    return (vector_index_ == rhs.getVectorIndex() &&
            mini_index_ == rhs.getMiniIndex());
  }
 
 private:
  NumSeqsType vector_index_ : 30;
  MiniIndex mini_index_ : 2;
};
/*--------------------------- NODE's INDEX -----------------------------------*/
 
/*--------------------------------------------------------------*/
/*--------------------------------------------------------------*/
const char REF_NAME[] = "REF";
const int MIN_NUM_TAXA = 3;
const RealNumType MIN_NEGATIVE =
    std::numeric_limits<RealNumType>::lowest();  // -FLT_MAX;
const RealNumType MIN_POSITIVE =
    (std::numeric_limits<RealNumType>::min)();  // FLT_MIN;
// const RealNumType INVERSE_MIN_POSITIVE = 1.0 / MIN_POSITIVE;
// const RealNumType LOG_MIN_POSITIVE = log(MIN_POSITIVE);
const RealNumType MAX_POSITIVE = (std::numeric_limits<RealNumType>::max)();
const RealNumType LOG_MAX_POSITIVE = log(MAX_POSITIVE);
 
template <class T>
constexpr T getMinCarryOver();
template <>
constexpr double getMinCarryOver<double>() {
  return 1e-250;
};  // of -308
template <>
constexpr float getMinCarryOver<float>() {
  return 1e-36f;
};  // of -38
 
const RealNumType MIN_CARRY_OVER = getMinCarryOver<RealNumType>();
const RealNumType MEAN_SUBS_PER_SITE = 0.02;
const RealNumType MAX_SUBS_PER_SITE = 0.067;
 
/*--------------------------------------------------------------*/
/*--------------------------------------------------------------*/
class Params {
 private:
  Params();
 
  // declare ParamsBuilder as a friend class to allow it access to the default
  // contructor Params();
  friend class ParamsBuilder;
 
 public:
  std::string aln_path;
 
  std::string ref_path;
 
  std::string ref_seqname;
 
  std::string aln_format_str;
 
  std::string sub_model_str;
 
  bool overwrite_output;
 
  bool fixed_blengths;
 
  RealNumType threshold_prob;
 
  RealNumType threshold_prob2;
 
  int failure_limit_sample;
 
  int failure_limit_subtree;
 
  int failure_limit_subtree_short_search;
 
  bool strict_stop_seeking_placement_sample;
 
  bool strict_stop_seeking_placement_subtree;
 
  bool strict_stop_seeking_placement_subtree_short_search;
 
  RealNumType thresh_log_lh_sample;
 
  RealNumType thresh_log_lh_subtree;
 
  RealNumType thresh_log_lh_subtree_short_search;
 
  RealNumType thresh_log_lh_failure;
 
  RealNumType fixed_min_blength;
 
  RealNumType min_blength_factor;
 
  RealNumType max_blength_factor;
 
  RealNumType min_blength_mid_factor;
 
  RealNumType thresh_diff_update;
 
  RealNumType thresh_diff_fold_update;
 
  PositionType mutation_update_period;
 
  std::string output_aln;
 
  std::string output_aln_format_str;
 
  std::string input_treefile;
 
  int32_t num_tree_improvement;
 
  RealNumType thresh_placement_cost;
 
  RealNumType thresh_entire_tree_improvement;
 
  RealNumType thresh_placement_cost_short_search;
 
  std::string tree_format_str;
 
  bool shallow_tree_search;
 
  char* output_testing;
 
  bool compute_aLRT_SH;
 
  PositionType aLRT_SH_replicates;
 
  RealNumType aLRT_SH_half_epsilon;
 
  uint32_t num_threads;
 
  uint64_t ran_seed;
 
  std::string output_prefix;
 
  bool allow_replace_input_tree;
 
  std::string seq_type_str;
 
  std::string tree_search_type_str;
 
  bool make_consistent;
    
  RealNumType max_subs_per_site;
    
  RealNumType mean_subs_per_site;
 
  /*
      TRUE to log debugging
   */
  // bool debug = false;
};
class ParamsBuilder {
 public:
  ParamsBuilder();
 
  ParamsBuilder& withRandomSeed(const uint64_t& seed);
 
  ParamsBuilder& withThreshProb(const double& thresh_prob);
 
  ParamsBuilder& withMinBlengthFactor(const double& min_blength_factor);
 
  ParamsBuilder& withMaxBlengthFactor(const double& max_blength_factor);
 
  ParamsBuilder& withFixedMinBlength(const double& fixed_min_blength);
 
  ParamsBuilder& withMutationUpdatePeriod(
      const int32_t& mutation_update_period);
 
  ParamsBuilder& withNumTreeTraversal(const int32_t& num_tree_traversal);
 
  ParamsBuilder& withSPRThresh(const double& SPR_thresh);
 
  ParamsBuilder& withStopTreeSearchThresh(const double& stop_search_thresh);
 
  std::unique_ptr<cmaple::Params> build();
 
 private:
  std::unique_ptr<Params> params_ptr{};
};
 
/*--------------------------------------------------------------*/
/*--------------------------------------------------------------*/
 
void outError(const char* error, bool quit = true);
 
void outError(const std::string& error, bool quit = true);
 
/*--------------------------------------------------------------*/
/*--------------------------------------------------------------*/
 
void outError(const char* error, const char* msg, bool quit = true);
 
void outError(const char* error, const std::string& msg, bool quit = true);
 
void outWarning(const char* warn);
void outWarning(const std::string& warn);
 
std::istream& safeGetline(std::istream& is, std::string& t);
 
inline bool controlchar(char& ch) {
  return ch <= 32;
}
 
inline bool is_newick_token(char& ch) {
  return ch == ':' || ch == ';' || ch == ',' || ch == ')' || ch == '(' ||
         ch == '[' || ch == ']';
}
 
bool renameString(std::string& name);
 
/*--------------------------------------------------------------*/
/*--------------------------------------------------------------*/
/*
        Error messages
 */
const char ERR_NO_TAXON[] = "Find no taxon with name ";
const char ERR_NO_AREA[] = "Find no area with name ";
const char ERR_NO_MEMORY[] = "Not enough memory!";
const char ERR_NEG_BRANCH[] = "Negative branch length not allowed.";
const char ERR_READ_INPUT[] =
    "File not found or incorrect input, pls check it again.";
const char ERR_READ_ANY[] =
    "Unidentified error while reading file, pls check it carefully again.";
 
/*--------------------------------------------------------------*/
/*--------------------------------------------------------------*/
 
std::string convertPosTypeToString(PositionType number);
std::string convertIntToString(int number);
std::string convertInt64ToString(int64_t number);
 
std::string convertDoubleToString(RealNumType number);
 
std::string convertDoubleToString(RealNumType number, uint8_t precision);
 
bool iEquals(const std::string& a, const std::string& b);
 
bool copyFile(const char SRC[], const char DEST[]);
 
bool fileExists(const std::string& strFilename);
 
int isDirectory(const char* path);
 
int convert_int(const char* str);
 
int64_t convert_int64(const char* str);
 
int convert_int(const char* str, int& end_pos);
 
void convert_int_vec(const char* str, IntVector& vec);
 
int64_t convert_int64(const char* str);
 
int64_t convert_int64(const char* str, int& end_pos);
 
RealNumType convert_real_number(const char* str);
 
RealNumType convert_real_number(const char* str, int& end_pos);
 
void convert_real_numbers(RealNumType*& arr, std::string input_str);
 
void convert_real_number_vec(const char* str,
                             RealNumberVector& vec,
                             char separator = ',');
 
void normalize_frequencies_from_index(RealNumType* freqs,
                                      int num_states,
                                      int starting_index);
 
inline void normalize_arr(RealNumType* const entries,
                          const int num_entries,
                          RealNumType sum_entries) {
  assert(num_entries > 0);
 
#ifndef NDEBUG
  // if (fabs(sum_entries) < 1e-5)
  //   outError("Sum of entries must be greater than zero!");
#endif
 
  sum_entries = 1.0 / sum_entries;
  for (int i = 0; i < num_entries; ++i)
    entries[i] *= sum_entries;
    // entries[i] /= sum_entries;
}
 
inline void normalize_arr(RealNumType* const entries, const int num_entries) {
  RealNumType sum_entries = 0;
  for (int i = 0; i < num_entries; ++i)
    sum_entries += entries[i];
  normalize_arr(entries, num_entries, sum_entries);
}
std::string convert_time(const RealNumType sec);
 
void convert_range(const char* str, int& lower, int& upper, int& step_size);
 
void convert_range(const char* str,
                   RealNumType& lower,
                   RealNumType& upper,
                   RealNumType& step_size);
 
void reinitDoubleArr(RealNumType*& arr,
                     StateType size,
                     bool delete_first = true,
                     bool set_zero = true);
 
void convert_string_vec(const char* str,
                        StrVector& str_vec,
                        char separator = ',');
 
PositionType convert_positiontype(const char* str);
 
bool is_number(const std::string& s);
 
void parseArg(int argc, char* argv[], Params& params);
 
void quickStartGuide();
 
void usage_cmaple();
 
void trimString(std::string& str);
 
int countPhysicalCPUCores();
 
template <class T1, class T2>
void quicksort(T1* arr, int left, int right, T2* arr2 = NULL) {
  if (left > right)
    return;
  assert(left <= right);
  int i = left, j = right;
  T1 pivot = arr[(left + right) / 2];
 
  /* partition */
  while (i <= j) {
    while (arr[i] < pivot)
      ++i;
    while (arr[j] > pivot)
      --j;
    if (i <= j) {
      T1 tmp = arr[i];
      arr[i] = arr[j];
      arr[j] = tmp;
      if (arr2) {
        T2 tmp2 = arr2[i];
        arr2[i] = arr2[j];
        arr2[j] = tmp2;
      }
      ++i;
      --j;
    }
  };
 
  /* recursion */
  if (left < j)
    quicksort(arr, left, j, arr2);
  if (i < right)
    quicksort(arr, i, right, arr2);
}
 
void usage_cmaple(char* argv[], bool full_command);
 
void printCopyright(std::ostream& out);
 
void setNumThreads(const int num_threads);
 
template <typename T, typename... Args>
std::unique_ptr<T> make_unique(Args&&... args) {
  return std::unique_ptr<T>(new T(std::forward<Args>(args)...));
}
 
void resetStream(std::istream& instream);
}  // namespace cmaple