utils.h

Go to the documentation of this file.

 
#include <armadillo>
 
#ifndef DISABLE_EIGEN
#include <Eigen/Core>
#endif 
 
#include <vector>
#include <random>
#include <algorithm>
#include <unordered_set>
#include <regex>
 
#include <boost/algorithm/string.hpp>
 
#ifndef UTILS_H_
#define UTILS_H_
 
namespace gds {

 
template<typename T>
static vector<double> GP_to_dosage( const vector<T> &v, const bool &missingToMean) {
  vector<double> res( v.size() / 3.0);
  vector<int> missing;
 
  // initialize
  int runningSum = 0, nValid = 0;
  double value;
 
  // for each entry in result
  // use two adjacent values
  for(int i=0; i<res.size(); i++){
    // compute dosage from genotype probabilties
    value = v[3*i]*0 + v[3*i+1]*1 + v[3*i+2]*2;
 
    // -9 is the missing value, so -18 is diploid
    if( value == -18){ 
      // if missing, set to NaN
      value = std::numeric_limits<double>::quiet_NaN();
      missing.push_back(i);
    }else{
      // for computing mean
      runningSum += value;
      nValid++;
    }
 
    // set dosage value
    res[i] = value;
  }
 
  // mean excluding NaNs
  double mu = runningSum / (double) nValid;
 
  // if missing values should be set to mean
  if( missingToMean ){
    // for each entry with a missing value, set to mean
    for(const int& i : missing) res[i] = mu;
  }
 
  return res;
}

static vector<double> intToDosage( const vector<int> &v, const bool &missingToMean) {
 
  // store and return result
  vector<double> res( v.size() / 2.0);
  vector<int> missing;
 
  // initialize
  int runningSum = 0, nValid = 0;
  double value;
 
  // for each entry in result
  // use two adjacent values
  for(int i=0; i<res.size(); i++){
    value = v[2*i] + v[2*i+1];
 
    // -9 is the missing value, so -18 is diploid
    if( value == -18){ 
      // if missing, set to NaN
      value = std::numeric_limits<double>::quiet_NaN();
      missing.push_back(i);
    }else{
      // for computing mean
      runningSum += value;
      nValid++;
    }
 
    // set dosage value
    res[i] = value;
  }
 
  // mean excluding NaNs
  double mu = runningSum / (double) nValid;
 
  // if missing values should be set to mean
  if( missingToMean ){
    // for each entry with a missing value, set to mean
    for(const int& i : missing) res[i] = mu;
  }
 
  return res;
}
 
template<typename T>
static size_t removeDuplicates(vector<T>& vec){
    unordered_set<T> seen;
 
    auto newEnd = remove_if(vec.begin(), vec.end(), [&seen](const T& value)
    {
        if (seen.find(value) != end(seen))
            return true;
 
        seen.insert(value);
        return false;
    });
 
    vec.erase(newEnd, vec.end());
 
    return vec.size();
}
 
 
static arma::vec colSums( const arma::mat &X){
 
  // row vector of 1's
  arma::rowvec ONE(X.n_rows, arma::fill::ones);
 
  // matrix multiplication to get sums
  arma::mat tmp = ONE * X;
 
  return arma::conv_to<arma::vec>::from( tmp );
}
 

static void standardize( arma::mat &X, const bool &center = true, const bool &scale = true, const double tol = 1e-10 ){
    
  double sqrt_rdf = sqrt(X.n_rows - 1.0);
 
  // if center, subtract mean of each column
  // if scale, divide by sd of each column
  // Note, norm() does not center the column
  //   this give results consistent with base::scale()
  //   when scale is FALSE
  for(size_t j=0; j<X.n_cols; j++){
    if( center ) X.col(j) -= mean(X.col(j));
    double sd = norm(X.col(j)) / sqrt_rdf;
    if( scale && sd > tol )  X.col(j) /= sd;
  }
}
 
#ifndef DISABLE_EIGEN
static void standardize( Eigen::MatrixXd &X, const bool &center = true, const bool &scale = true, const double tol = 1e-10 ){
    
  double sqrt_rdf = sqrt(X.rows() - 1.0);
  if(center) X.rowwise() -= X.colwise().mean(); // centering
  if(scale) {
    // if X is centered, then we can convert norm to sd
    for(size_t j=0; j<X.cols(); j++) {
      double sd = X.col(j).norm() / sqrt_rdf ; // sd
      if(sd > tol) X.col(j) /= sd;
    }
  }
}
 
 
static void standardize_rows( Eigen::MatrixXd &X, const bool &center = true, const bool &scale = true, const double tol = 1e-10 ){
    
  double sqrt_rdf = sqrt(X.cols() - 1.0);
  if(center) X.colwise() -= X.rowwise().mean(); // centering
  if(scale) {
    // if X is centered, then we can convert norm to sd
    for(size_t j=0; j<X.rows(); j++) {
      double sd = X.row(j).norm() / sqrt_rdf ; // sd
      if(sd > tol) X.row(j) /= sd;
    }
  }
}
#endif

static bool isOnlyDigits(const std::string& s){
  int n = count_if(s.begin(), s.end(),
                       [](unsigned char c){ return isdigit(c); } 
                      );
  return( n == s.size());
}
 
 
static void nanToMean( arma::vec & v){
    // get indeces of finite elements
  arma::uvec idx = arma::find_finite(v);
 
  // if number of finite elements is less than the total
  if( idx.n_elem < v.n_elem ){
    // compute mean from finite elements
    double mu = arma::mean( v.elem(idx));
 
    // replace nan with mu
    v.replace(arma::datum::nan, mu);
  }
}
 
 
typedef enum {
  VCF,
  VCFGZ,
  BCF,
  BGEN,
  PGEN,
  PBED,
  OTHER
} FileType;
 
static string toString( FileType x){
 
  switch(x){
    case VCF:   return "vcf";
    case VCFGZ: return "vcf.gz";
    case BCF:   return "bcf";
    case BGEN:  return "bgen";
    case PGEN:  return "pgen";
    case PBED:  return "bed";
    case OTHER: return "other";
    default:    return "other";
  }
}
 
 
 
static FileType getFileType( const string &file ){
 
  FileType ft = OTHER;
 
  if( regex_search( file, regex("\\.vcf$")) ){
    ft = VCF;
  }else if( regex_search( file, regex("\\.vcf\\.gz$")) ){
    ft = VCFGZ;
  }else if( regex_search( file, regex("\\.bcf$")) ) {
    ft = BCF;
  }else if( regex_search( file, regex("\\.bgen$")) ){
    ft = BGEN;
  }else if( regex_search( file, regex("\\.pgen$")) ){
    ft = PGEN;
  } if( regex_search( file, regex("\\.bed$")) ){
    ft = PBED;
  }
 
  return ft;
}
 
 
 
template<typename T>
static vector<T> cast_elements( const vector<string> &v ){
  vector<T> output(0, v.size());
 
  for (auto &s : v) {
    stringstream parser(s);
    T x = 0;
    parser >> x;
    output.push_back(x);
  }
  return output;
}
 
 
static vector<int> stoi_vec( const vector<string> &v ){
 
  vector<int> output;
  output.reserve(v.size());
 
  for (const string& s : v) {
    output.push_back(std::stoi(s));
  }
  return output;
}
 
 

static vector<string> splitRegionString( string regionString){
 
  vector<string> regions;
 
  // regionString is string of chr:start-end delim by "\t,\n"
  // remove spaces, then split based on delim
  boost::erase_all(regionString, " ");
  boost::split(regions, regionString, boost::is_any_of("\t,\n"));
 
  // remove duplicate regions, but preserve order
  removeDuplicates( regions );
 
  return regions;
}
 

template<typename T, typename T2>
static vector<T> subset_vector(const vector<T> &x, const vector<T2> &idx){
 
  // initialize x_subset to have size idx.size()
  vector<T> x_subset;
  x_subset.reserve(idx.size());
 
  for(auto i: idx){
      if( i > x.size() ){
          throw std::out_of_range("Index is out of bounds"); 
      }
      x_subset.push_back( x[i] );
  }
 
  return x_subset;
}
 
 
template<typename T>
static void print_vec(const string &title, const vector<T> &x){
  Rcpp::Rcout << title << "\n";
  for(auto a: x) Rcpp::Rcout << a << " ";
  Rcpp::Rcout << endl;
}
 
 
template<typename T>
std::vector<std::vector<T>> chunk_vector(const std::vector<T>& vec, int k) {
  std::vector<std::vector<T> > result;
  int size = vec.size();
  int chunk_size = size / k;
  int remainder = size % k;
 
  int start = 0;
  for (int i = 0; i < k; ++i) {
    int current_chunk_size = chunk_size;
    if (i < remainder) {
        current_chunk_size++;
    }
    std::vector<T> chunk(vec.begin() + start, vec.begin() + start + current_chunk_size);
    result.push_back(chunk);
    start += current_chunk_size;
  }
  return result;
}
 
template<typename T>
static vector<unsigned int> which_in( const vector<T> &v1, const vector<T> &v2){
 
  // Create a hash map (unordered_map) of elements in v2 for efficient lookups
  unordered_map<T, bool> v2_elements;
  for (auto &x : v2) {
      v2_elements[x] = true;
  }
 
  // Iterate through v1 and check if the element exists in the hash map   
  vector<unsigned int> shared_indices;
  for (int i = 0; i < v1.size(); ++i) {
      // Check if the element exists in the map
      if (v2_elements.count(v1[i])) { 
        // Add the index from v1
          shared_indices.push_back(i); 
      }
  }
 
  return shared_indices;
}
 
 
 
 
 
 
} // end namespace
#endif