math.h

Go to the documentation of this file.

 
#pragma once
 
#include "jit_normalization.h"
#include "threads.h"
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <numeric>
#include <vector>
 
namespace slope {
 
template<typename T>
int
sign(T val)
{
  return (T(0) < val) - (val < T(0));
}
 
template<typename T>
Eigen::ArrayXd
cumSum(const T& x)
{
  std::vector<double> cum_sum(x.size());
  std::partial_sum(x.begin(), x.end(), cum_sum.begin(), std::plus<double>());
 
  Eigen::Map<Eigen::ArrayXd> out(cum_sum.data(), cum_sum.size());
 
  return out;
}
 
template<typename T>
T
sigmoid(const T& x)
{
  return 1.0 / (1.0 + std::exp(-x));
}
 
template<typename T>
T
logit(const T& x)
{
  assert(x > 0 && x < 1 && "Input must be in (0, 1)");
 
  return std::log(x) - std::log1p(-x);
}
 
template<typename T>
T
clamp(const T& x, const T& lo, const T& hi)
{
  return x < lo ? lo : x > hi ? hi : x;
}
 
Eigen::VectorXd
logSumExp(const Eigen::MatrixXd& a);
 
Eigen::MatrixXd
softmax(const Eigen::MatrixXd& x);
 
template<typename T>
Eigen::MatrixXd
linearPredictor(const T& x,
                const std::vector<int>& active_set,
                const Eigen::VectorXd& beta0,
                const Eigen::VectorXd& beta,
                const Eigen::VectorXd& x_centers,
                const Eigen::VectorXd& x_scales,
                const JitNormalization jit_normalization,
                const bool intercept)
{
  int n = x.rows();
  int p = x.cols();
  int m = beta0.size();
 
  Eigen::MatrixXd eta = Eigen::MatrixXd::Zero(n, m);
 
#ifdef _OPENMP
  bool large_problem = active_set.size() > 100 && n * active_set.size() > 1e7;
#pragma omp parallel num_threads(Threads::get()) if (large_problem)
#endif
  {
    Eigen::MatrixXd eta_local = Eigen::MatrixXd::Zero(n, m);
 
#ifdef _OPENMP
#pragma omp for nowait
#endif
    for (size_t i = 0; i < active_set.size(); ++i) {
      int ind = active_set[i];
      auto [k, j] = std::div(ind, p);
 
      switch (jit_normalization) {
        case JitNormalization::Both:
          eta_local.col(k) += x.col(j) * beta(ind) / x_scales(j);
          eta_local.col(k).array() -= beta(ind) * x_centers(j) / x_scales(j);
          break;
 
        case JitNormalization::Center:
          eta_local.col(k) += x.col(j) * beta(ind);
          eta_local.col(k).array() -= beta(ind) * x_centers(j);
          break;
 
        case JitNormalization::Scale:
          eta_local.col(k) += x.col(j) * beta(ind) / x_scales(j);
          break;
 
        case JitNormalization::None:
          eta_local.col(k) += x.col(j) * beta(ind);
          break;
      }
    }
 
#ifdef _OPENMP
#pragma omp critical
#endif
    {
      eta += eta_local;
    }
  }
 
  if (intercept) {
    eta.rowwise() += beta0.transpose();
  }
 
  return eta;
}
 
template<typename T>
void
updateGradient(Eigen::VectorXd& gradient,
               const T& x,
               const Eigen::MatrixXd& residual,
               const std::vector<int>& active_set,
               const Eigen::VectorXd& x_centers,
               const Eigen::VectorXd& x_scales,
               const Eigen::VectorXd& w,
               const JitNormalization jit_normalization)
{
  const int n = x.rows();
  const int p = x.cols();
  const int m = residual.cols();
 
  assert(gradient.size() == p * m &&
         "Gradient matrix has incorrect dimensions");
 
  Eigen::MatrixXd weighted_residual(n, m);
  Eigen::ArrayXd wr_sums(m);
 
#ifdef _OPENMP
  bool large_problem = active_set.size() > 100 && n * active_set.size() > 1e5;
#pragma omp parallel for num_threads(Threads::get()) if (large_problem)
#endif
  for (int k = 0; k < m; ++k) {
    weighted_residual.col(k) = residual.col(k).cwiseProduct(w);
    wr_sums(k) = weighted_residual.col(k).sum();
  }
 
#ifdef _OPENMP
#pragma omp parallel for num_threads(Threads::get()) if (large_problem)
#endif
  for (size_t i = 0; i < active_set.size(); ++i) {
    int ind = active_set[i];
    auto [k, j] = std::div(ind, p);
 
    switch (jit_normalization) {
      case JitNormalization::Both:
        gradient(ind) =
          (x.col(j).dot(weighted_residual.col(k)) - x_centers(j) * wr_sums(k)) /
          (x_scales(j) * n);
        break;
      case JitNormalization::Center:
        gradient(ind) =
          (x.col(j).dot(weighted_residual.col(k)) - x_centers(j) * wr_sums(k)) /
          n;
        break;
      case JitNormalization::Scale:
        gradient(ind) =
          x.col(j).dot(weighted_residual.col(k)) / (x_scales(j) * n);
        break;
      case JitNormalization::None:
        gradient(ind) = x.col(j).dot(weighted_residual.col(k)) / n;
        break;
    }
  }
}
 
template<typename T>
void
offsetGradient(Eigen::VectorXd& gradient,
               const T& x,
               const Eigen::VectorXd& offset,
               const std::vector<int>& active_set,
               const Eigen::VectorXd& x_centers,
               const Eigen::VectorXd& x_scales,
               const JitNormalization jit_normalization)
{
  const int n = x.rows();
  const int p = x.cols();
 
  for (size_t i = 0; i < active_set.size(); ++i) {
    int ind = active_set[i];
    auto [k, j] = std::div(ind, p);
 
    switch (jit_normalization) {
      case JitNormalization::Both:
        gradient(ind) -=
          offset(k) * (x.col(j).sum() / n - x_centers(j)) / x_scales(j);
        break;
      case JitNormalization::Center:
        gradient(ind) -= offset(k) * (x.col(j).sum() / n - x_centers(j));
        break;
      case JitNormalization::Scale:
        gradient(ind) -= offset(k) * x.col(j).sum() / (n * x_scales(j));
        break;
      case JitNormalization::None:
        gradient(ind) -= offset(k) * x.col(j).sum() / n;
        break;
    }
  }
}
 
std::vector<int>
setUnion(const std::vector<int>& a, const std::vector<int>& b);
 
std::vector<int>
setDiff(const std::vector<int>& a, const std::vector<int>& b);
 
template<typename T>
int
whichMax(const T& x)
{
  return std::distance(x.begin(), std::max_element(x.begin(), x.end()));
}
 
template<typename T>
int
whichMin(const T& x)
{
  return std::distance(x.begin(), std::min_element(x.begin(), x.end()));
}
 
template<typename T, typename Comparator>
int
whichBest(const T& x, const Comparator& comp)
{
  return std::distance(x.begin(), std::max_element(x.begin(), x.end(), comp));
}
 
Eigen::ArrayXd
geomSpace(const double start, const double end, const int n);
 
template<typename T>
Eigen::VectorXd
l1Norms(const T& x)
{
  const int p = x.cols();
 
  Eigen::VectorXd out(p);
 
  for (int j = 0; j < p; ++j) {
    out(j) = x.col(j).cwiseAbs().sum();
  }
 
  return out;
}
 
Eigen::VectorXd
l2Norms(const Eigen::SparseMatrix<double>& x);
 
Eigen::VectorXd
l2Norms(const Eigen::MatrixXd& x);
 
Eigen::VectorXd
maxAbs(const Eigen::SparseMatrix<double>& x);
 
Eigen::VectorXd
maxAbs(const Eigen::MatrixXd& x);
 
Eigen::VectorXd
means(const Eigen::SparseMatrix<double>& x);
 
Eigen::VectorXd
means(const Eigen::MatrixXd& x);
 
Eigen::VectorXd
stdDevs(const Eigen::SparseMatrix<double>& x);
 
Eigen::VectorXd
stdDevs(const Eigen::MatrixXd& x);
 
Eigen::VectorXd
ranges(const Eigen::SparseMatrix<double>& x);
 
Eigen::VectorXd
ranges(const Eigen::MatrixXd& x);
 
Eigen::VectorXd
mins(const Eigen::SparseMatrix<double>& x);
 
Eigen::VectorXd
mins(const Eigen::MatrixXd& x);
 
} // namespace slope