math.h

Go to the documentation of this file.

 
#pragma once
 
#include "clusters.h"
#include "jit_normalization.h"
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <numeric>
#include <vector>
 
#ifdef _OPENMP
#include "threads.h"
#endif
 
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 (int i = 0; i < static_cast<int>(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 (int i = 0; i < static_cast<int>(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);
 
template<typename T>
Eigen::VectorXd
clusterGradient(Eigen::VectorXd& beta,
                Eigen::VectorXd& residual,
                Clusters& clusters,
                const T& x,
                const Eigen::VectorXd& w,
                const Eigen::VectorXd& x_centers,
                const Eigen::VectorXd& x_scales,
                const JitNormalization jit_normalization)
{
  using namespace Eigen;
 
  const int n = x.rows();
  const int n_clusters = clusters.n_clusters();
 
  Eigen::VectorXd gradient = Eigen::VectorXd::Zero(n_clusters);
 
  for (int j = 0; j < n_clusters; ++j) {
    double c_old = clusters.coeff(j);
 
    if (c_old == 0) {
      gradient(j) = 0;
      continue;
    }
 
    int cluster_size = clusters.cluster_size(j);
    std::vector<int> s;
    s.reserve(cluster_size);
 
    for (auto c_it = clusters.cbegin(j); c_it != clusters.cend(j); ++c_it) {
      double s_k = sign(beta(*c_it));
      s.emplace_back(s_k);
    }
 
    double hessian_j = 1;
    double gradient_j = 0;
 
    if (cluster_size == 1) {
      int k = *clusters.cbegin(j);
      std::tie(gradient_j, hessian_j) = computeGradientAndHessian(
        x, k, w, residual, x_centers, x_scales, s[0], jit_normalization, n);
    } else {
      std::tie(hessian_j, gradient_j) = computeClusterGradientAndHessian(
        x, j, s, clusters, w, residual, x_centers, x_scales, jit_normalization);
    }
 
    gradient(j) = gradient_j;
  }
 
  return gradient;
}
 
} // namespace slope