feat: add initial c++ source for libslope

src/main.cpp

@@ -0,0 +1,49 @@
+#include "slope/slope.h"
+#include <Eigen/Core>
+#include <Eigen/LU>
+#include <pybind11/eigen.h>
+#include <pybind11/pybind11.h>
+
+#define STRINGIFY(x) #x
+#define MACRO_STRINGIFY(x) STRINGIFY(x)
+
+int
+fit_slope(const Eigen::MatrixXd x,
+          const Eigen::MatrixXd y,
+          Eigen::ArrayXd lambda,
+          Eigen::ArrayXd alpha)
+{
+  auto result = slope::slope(x, y, alpha, lambda);
+
+  return 0;
+}
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(_sortedl1, m)
+{
+  m.doc() = R"pbdoc(
+        Pybind11 example plugin
+        -----------------------
+
+        .. currentmodule:: _sortedl1
+
+        .. autosummary::
+           :toctree: _generate
+
+           fit_slope
+    )pbdoc";
+
+  m.def("fit_slope", &fit_slope, R"pbdoc(
+        Fit SLOPE
+        ---------------
+
+        Let's fit SLOPE!
+    )pbdoc");
+
+#ifdef VERSION_INFO
+  m.attr("__version__") = MACRO_STRINGIFY(VERSION_INFO);
+#else
+  m.attr("__version__") = "dev";
+#endif
+}

src/slope/cd.h

@@ -0,0 +1,138 @@
+#pragma once
+
+#include "clusters.h"
+#include "slope_threshold.h"
+#include "sorted_l1_norm.h"
+#include <Eigen/Core>
+#include <vector>
+
+namespace slope {
+
+template<typename T>
+void
+coordinateDescent(double& beta0,
+                  Eigen::VectorXd& beta,
+                  Eigen::VectorXd& residual,
+                  Clusters& clusters,
+                  const T& x,
+                  const Eigen::VectorXd& w,
+                  const Eigen::VectorXd& z,
+                  const SortedL1Norm& sl1_norm,
+                  const Eigen::VectorXd& x_centers,
+                  const Eigen::VectorXd& x_scales,
+                  bool intercept,
+                  bool standardize,
+                  bool update_clusters,
+                  int print_level)
+{
+  using namespace Eigen;
+
+  const int n = x.rows();
+
+  for (int j = 0; j < clusters.n_clusters(); ++j) {
+    double c_old = clusters.coeff(j);
+
+    if (c_old == 0) {
+      // We do not update the zero cluster because it can be very large, but
+      // often does not change.
+      continue;
+    }
+
+    std::vector<int> s;
+    int cluster_size = clusters.cluster_size(j);
+    s.reserve(cluster_size);
+
+    double hessian_j = 1;
+    double gradient_j = 0;
+    VectorXd x_s(n);
+
+    if (cluster_size == 1) {
+      int k = *clusters.cbegin(j);
+      double s_k = sign(beta(k));
+      s.emplace_back(s_k);
+
+      if (standardize) {
+        gradient_j = -s_k *
+                     (x.col(k).cwiseProduct(w).dot(residual) -
+                      w.dot(residual) * x_centers(k)) /
+                     (n * x_scales(k));
+        // TODO: Consider caching the hessian values. We need to invalidate the
+        // cache every time the cluster is updated or the  signs flip relatively
+        // inside the cluster.
+        hessian_j =
+          (x.col(k).cwiseAbs2().dot(w) - 2 * x_centers(k) * x.col(k).dot(w) +
+           std::pow(x_centers(k), 2) * w.sum()) /
+          (std::pow(x_scales(k), 2) * n);
+      } else {
+        gradient_j = -s_k * (x.col(k).cwiseProduct(w).dot(residual)) / n;
+        hessian_j = x.col(k).cwiseAbs2().dot(w) / n;
+      }
+    } else {
+      // There's no reasonable just-in-time standardization approach for sparse
+      // design matrices when there are clusters in the data, so we need to
+      // reduce to a dense column vector.
+      x_s.setZero();
+
+      for (auto c_it = clusters.cbegin(j); c_it != clusters.cend(j); ++c_it) {
+        int k = *c_it;
+        double s_k = sign(beta(k));
+        s.emplace_back(s_k);
+
+        if (standardize) {
+          x_s += x.col(k) * (s_k / x_scales(k));
+          x_s.array() -= x_centers(k) * s_k / x_scales(k);
+        } else {
+          x_s += x.col(k) * s_k;
+        }
+      }
+
+      hessian_j = x_s.cwiseAbs2().dot(w) / n;
+      gradient_j = -x_s.cwiseProduct(w).dot(residual) / n;
+    }
+
+    auto thresholding_results =
+      slopeThreshold(c_old - gradient_j / hessian_j,
+                     j,
+                     sl1_norm.lambda * sl1_norm.getAlpha() / hessian_j,
+                     clusters);
+
+    double c_tilde = thresholding_results.value;
+    int new_index = thresholding_results.new_index;
+
+    auto s_it = s.cbegin();
+    auto c_it = clusters.cbegin(j);
+    for (; c_it != clusters.cend(j); ++c_it, ++s_it) {
+      beta(*c_it) = c_tilde * (*s_it);
+    }
+
+    double c_diff = c_old - c_tilde;
+
+    if (c_diff != 0) {
+      if (cluster_size == 1) {
+        int k = *clusters.cbegin(j);
+
+        if (standardize) {
+          residual += x.col(k) * (s[0] * c_diff / x_scales(k));
+          residual.array() -= x_centers(k) * s[0] * c_diff / x_scales(k);
+        } else {
+          residual += x.col(k) * (s[0] * c_diff);
+        }
+      } else {
+        residual += x_s * c_diff;
+      }
+    }
+
+    if (update_clusters) {
+      clusters.update(j, new_index, std::abs(c_tilde));
+    } else {
+      clusters.setCoeff(j, std::abs(c_tilde));
+    }
+
+    if (intercept) {
+      double beta0_update = residual.dot(w) / w.sum();
+      residual.array() -= beta0_update;
+      beta0 += beta0_update;
+    }
+  }
+}
+}

src/slope/clusters.cpp

@@ -0,0 +1,206 @@
+#include "clusters.h"
+
+namespace slope {
+
+Clusters::Clusters(const Eigen::VectorXd& beta)
+{
+  update(beta);
+}
+
+std::vector<int>::iterator
+Clusters::begin(const int i)
+{
+  return c_ind.begin() + c_ptr[i];
+}
+
+std::vector<int>::iterator
+Clusters::end(const int i)
+{
+  return c_ind.begin() + c_ptr[i + 1];
+}
+
+std::vector<int>::const_iterator
+Clusters::cbegin(const int i) const
+{
+  return c_ind.cbegin() + c_ptr[i];
+}
+
+std::vector<int>::const_iterator
+Clusters::cend(const int i) const
+{
+  return c_ind.cbegin() + c_ptr[i + 1];
+}
+
+int
+Clusters::cluster_size(const int i) const
+{
+  return c_ptr[i + 1] - c_ptr[i];
+}
+
+int
+Clusters::pointer(const int i) const
+{
+  return c_ptr[i];
+}
+
+int
+Clusters::n_clusters() const
+{
+  return c.size();
+}
+
+double
+Clusters::coeff(const int i) const
+{
+  return c[i];
+}
+
+void
+Clusters::setCoeff(const int i, const double x)
+{
+  c[i] = x;
+}
+
+std::vector<double>
+Clusters::coeffs() const
+{
+  return c;
+}
+
+std::vector<int>
+Clusters::indices() const
+{
+  return c_ind;
+}
+
+std::vector<int>
+Clusters::pointers() const
+{
+  return c_ptr;
+}
+
+void
+Clusters::update(const int old_index, const int new_index, const double c_new)
+{
+  auto c_old = coeff(old_index);
+
+  if (c_new != c_old) {
+    if (c_new == coeff(new_index)) {
+      merge(old_index, new_index);
+    } else {
+      setCoeff(old_index, c_new);
+      if (old_index != new_index) {
+        reorder(old_index, new_index);
+      }
+    }
+  }
+}
+
+void
+Clusters::update(const Eigen::VectorXd& beta)
+{
+  using sort_pair = std::pair<double, int>;
+
+  c.clear();
+  c_ind.clear();
+  c_ptr.clear();
+
+  std::vector<sort_pair> sorted;
+  sorted.reserve(beta.size());
+
+  for (int i = 0; i < beta.size(); ++i) {
+    sorted.emplace_back(std::abs(beta(i)), i);
+  }
+
+  std::sort(sorted.begin(), sorted.end(), std::greater<sort_pair>());
+
+  c_ind.reserve(sorted.size());
+  for (const auto& sorted_i : sorted) {
+    c_ind.emplace_back(sorted_i.second);
+  }
+
+  std::vector<sort_pair> sorted_unique;
+  std::unique_copy(sorted.begin(),
+                   sorted.end(),
+                   std::back_inserter(sorted_unique),
+                   [](const sort_pair& a, const sort_pair& b) -> bool {
+                     return a.first == b.first;
+                   });
+
+  c.reserve(sorted_unique.size());
+  for (const auto& sorted_unique_i : sorted_unique) {
+    c.emplace_back(std::move(sorted_unique_i.first));
+  }
+
+  c_ptr.reserve(c.size() + 1);
+  c_ptr.emplace_back(0);
+
+  auto range_start = sorted.begin();
+  for (const auto& c_i : c) {
+    auto range_end = std::find_if(
+      std::next(range_start), sorted.end(), [&c_i](const sort_pair& x) -> bool {
+        return x.first != c_i;
+      });
+    c_ptr.emplace_back(std::distance(range_start, range_end));
+    range_start = range_end;
+  }
+
+  std::partial_sum(c_ptr.cbegin(), c_ptr.cend(), c_ptr.begin());
+}
+
+void
+Clusters::reorder(const int old_index, const int new_index)
+{
+  auto c_size = cluster_size(old_index);
+
+  // update coefficients
+  move_elements(c, old_index, new_index, 1);
+
+  // update indices
+  move_elements(c_ind, pointer(old_index), pointer(new_index), c_size);
+
+  // update pointers
+  if (new_index < old_index) {
+    move_elements(c_ptr, old_index + 1, new_index + 1, 1);
+
+    std::for_each(c_ptr.begin() + new_index + 1,
+                  c_ptr.begin() + old_index + 2,
+                  [c_size](int& x) { x += c_size; });
+
+    c_ptr[new_index + 1] = c_ptr[new_index] + c_size;
+  } else {
+    move_elements(c_ptr, old_index, new_index, 1);
+
+    std::for_each(c_ptr.begin() + old_index,
+                  c_ptr.begin() + new_index,
+                  [c_size](int& x) { x -= c_size; });
+    c_ptr[new_index] = c_ptr[new_index + 1] - c_size;
+  }
+}
+
+void
+Clusters::merge(const int old_index, const int new_index)
+{
+  auto c_size = cluster_size(old_index);
+
+  // update coefficients
+  c.erase(c.cbegin() + old_index);
+
+  // update indices
+  move_elements(c_ind, pointer(old_index), pointer(new_index), c_size);
+
+  // update pointers
+  if (new_index < old_index) {
+    std::for_each(c_ptr.begin() + new_index + 1,
+                  c_ptr.begin() + old_index + 1,
+                  [c_size](int& x) { x += c_size; });
+  } else {
+    std::for_each(c_ptr.begin() + old_index + 1,
+                  c_ptr.begin() + new_index + 1,
+                  [c_size](int& x) { x -= c_size; });
+  }
+
+  c_ptr.erase(c_ptr.begin() + old_index + 1);
+}
+
+} // namespace slope