utils.h

Go to the documentation of this file.

 
#pragma once
 
#include "eigen_compat.h"
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <algorithm>
#include <numeric>
#include <set>
#include <string>
#include <unordered_set>
#include <vector>
 
namespace slope {
 
using slope::all;
 
template<typename Derived>
Eigen::Index
nonZeros(const Eigen::SparseMatrixBase<Derived>& x)
{
  return x.derived().nonZeros();
}
 
template<typename Derived>
Eigen::Index
nonZeros(const Eigen::DenseBase<Derived>& x)
{
  return (x.derived().array() != 0).count();
}
 
template<typename T>
void
sort(T& v, const bool descending = false)
{
  if (descending) {
    std::sort(
      v.data(), v.data() + v.size(), std::greater<typename T::value_type>());
  } else {
    std::sort(
      v.data(), v.data() + v.size(), std::less<typename T::value_type>());
  }
}
 
template<typename T>
std::vector<int>
which(const T& x)
{
  std::vector<int> out;
  for (int i = 0; i < x.size(); i++) {
    if (x[i]) {
      out.emplace_back(i);
    }
  }
 
  return out;
}
 
template<typename T>
std::vector<int>
sortIndex(T& v, const bool descending = false)
{
  using namespace std;
 
  vector<int> idx(v.size());
  iota(idx.begin(), idx.end(), 0);
 
  if (descending) {
    sort(idx.begin(), idx.end(), [&v](int i, int j) { return v[i] > v[j]; });
  } else {
    sort(idx.begin(), idx.end(), [&v](int i, int j) { return v[i] < v[j]; });
  }
 
  return idx;
}
 
template<typename T>
void
permute(T& values, const std::vector<int>& ind)
{
  T out(values.size());
 
  for (int i = 0; i < values.size(); ++i)
    out[i] = std::move(values[ind[i]]);
 
  values = std::move(out);
}
 
template<typename T>
void
inversePermute(T& values, const std::vector<int>& ind)
{
  T out(values.size()); 
  for (int i = 0; i < values.size(); ++i)
    out[ind[i]] = std::move(values[i]);
 
  values = std::move(out);
}
 
template<typename T>
void
move_elements(std::vector<T>& v, const int from, const int to, const int size)
{
  assert(from >= 0);
  assert(to >= 0);
  assert(size >= 0);
  assert(from != to);
 
  if (from > to) {
    assert(from + size <= static_cast<int>(v.size()));
    std::rotate(v.begin() + to, v.begin() + from, v.begin() + from + size);
  } else {
    assert(to + size <= static_cast<int>(v.size()));
    std::rotate(
      v.begin() + from, v.begin() + from + size, v.begin() + to + size);
  }
}
 
void
validateOption(const std::string& value,
               const std::set<std::string>& valid_options,
               const std::string& parameter_name);
 
template<typename T>
T
subset(const Eigen::EigenBase<T>& x, const std::vector<int>& indices)
{
  return subset(x.derived(), indices);
}
 
template<typename T>
typename Eigen::MatrixBase<T>::PlainObject
subset(const Eigen::DenseBase<T>& x, const std::vector<int>& indices)
{
  return x.derived()(indices, all);
}
 
template<typename T>
T
subset(const Eigen::SparseMatrixBase<T>& x, const std::vector<int>& indices)
{
  std::vector<Eigen::Triplet<double>> triplets;
  triplets.reserve(slope::nonZeros(x.derived()));
 
  for (int j = 0; j < x.cols(); ++j) {
    for (typename T::InnerIterator it(x.derived(), j); it; ++it) {
      auto it_idx = std::find(indices.begin(), indices.end(), it.row());
 
      if (it_idx != indices.end()) {
        int new_row = std::distance(indices.begin(), it_idx);
        triplets.emplace_back(new_row, j, it.value());
      }
    }
  }
 
  T out(indices.size(), x.cols());
  out.setFromTriplets(triplets.begin(), triplets.end());
 
  return out;
}
 
template<typename T>
T
subsetCols(const Eigen::MatrixBase<T>& x, const std::vector<int>& indices)
{
  return x.derived()(all, indices);
}
 
template<typename T>
T
subsetCols(const Eigen::SparseMatrixBase<T>& x, const std::vector<int>& indices)
{
  std::vector<Eigen::Triplet<double>> triplets;
  triplets.reserve(slope::nonZeros(x.derived()));
 
  for (size_t j_idx = 0; j_idx < indices.size(); ++j_idx) {
    int j = indices[j_idx];
    for (typename T::InnerIterator it(x.derived(), j); it; ++it) {
      triplets.emplace_back(it.row(), j_idx, it.value());
    }
  }
 
  T out(x.rows(), indices.size());
  out.setFromTriplets(triplets.begin(), triplets.end());
 
  return out;
}
 
inline std::unordered_set<double>
unique(const Eigen::MatrixXd& x)
{
  std::unordered_set<double> unique;
  for (Eigen::Index j = 0; j < x.cols(); j++) {
    for (Eigen::Index i = 0; i < x.rows(); i++) {
      unique.insert(x(i, j));
    }
  }
 
  return unique;
}
 
} // namespace slope