cv.h

Go to the documentation of this file.

 
#pragma once
 
#include "eigen_compat.h"
#include "folds.h"
#include "score.h"
#include "slope.h"
#include <atomic>
#include <vector>
 
#ifdef _OPENMP
#include <omp.h>
#endif
 
namespace slope {
 
using slope::all;
 
struct GridResult
{
  Eigen::MatrixXd score;
 
  std::map<std::string, double> params;
 
  Eigen::ArrayXd alphas;
 
  Eigen::ArrayXd mean_scores;
 
  Eigen::ArrayXd std_errors;
};
 
struct CvResult
{
  std::vector<GridResult> results;
 
  std::map<std::string, double> best_params;
 
  double best_score;
 
  int best_ind;
 
  int best_alpha_ind;
};
 
struct CvConfig
{
  int n_folds = 10;
 
  int n_repeats = 1;
 
  std::string metric = "mse";
 
  uint64_t random_seed = 42;
 
  bool copy_x = true;
 
  std::map<std::string, std::vector<double>> hyperparams;
 
  std::map<std::string, std::vector<double>> default_hyperparams = {
    { "q", { 0.1 } },
    { "gamma", { 0.0 } },
  };
 
  std::optional<std::vector<std::vector<std::vector<int>>>> predefined_folds;
};
 
namespace detail {
 
std::vector<std::map<std::string, double>>
createGrid(const std::map<std::string, std::vector<double>>& param_values);
 
void
findBestParameters(CvResult& cv_result, const std::unique_ptr<Score>& scorer);
 
template<typename T>
Eigen::ArrayXd
fitToFold(Eigen::MatrixBase<T>& x,
          const Eigen::MatrixXd& y,
          const Folds& folds,
          const std::unique_ptr<Loss>& loss,
          const std::unique_ptr<Score>& scorer,
          const Eigen::ArrayXd& alphas,
          Slope& thread_model,
          const int fold,
          const int rep,
          const double gamma = 0.0,
          const bool copy_x = true)
{
  Eigen::ArrayXd scores = Eigen::ArrayXd::Zero(alphas.size());
 
  if (copy_x) {
    thread_model.setModifyX(true);
 
    auto [x_train, y_train, x_test, y_test] = folds.split(x, y, fold, rep);
 
    auto path = thread_model.path(x_train, y_train, alphas);
 
    if (gamma > 0) {
      path = thread_model.relax(path, x_train, y_train, gamma);
    }
 
    for (int j = 0; j < path.size(); ++j) {
      auto eta = path(j).predict(x_test, "linear");
      scores(j) = scorer->eval(eta, y_test, loss);
    }
 
  } else {
    thread_model.setModifyX(false);
 
    auto train_idx = folds.getTrainingIndices(fold, rep);
    auto test_idx = folds.getTestIndices(fold, rep);
 
    // Create views
    auto x_train = x(train_idx, all);
    auto x_test = x(test_idx, all);
 
    Eigen::MatrixXd y_train = y(train_idx, all);
    Eigen::MatrixXd y_test = y(test_idx, all);
 
    auto path = thread_model.path(x_train, y_train, alphas);
 
    if (gamma > 0) {
      path = thread_model.relax(path, x_train, y_train, gamma);
    }
 
    for (int j = 0; j < path.size(); ++j) {
      auto eta = path(j).predict(x_test, "linear");
      scores(j) = scorer->eval(eta, y_test, loss);
    }
  }
 
  return scores;
}
 
template<typename T>
Eigen::ArrayXd
fitToFold(Eigen::SparseMatrixBase<T>& x,
          const Eigen::MatrixXd& y,
          const Folds& folds,
          const std::unique_ptr<Loss>& loss,
          const std::unique_ptr<Score>& scorer,
          const Eigen::ArrayXd& alphas,
          Slope& thread_model,
          const int fold,
          const int rep,
          const double gamma = 0.0,
          const bool copy_x = true)
{
  thread_model.setModifyX(true);
 
  auto [x_train, y_train, x_test, y_test] = folds.split(x, y, fold, rep);
 
  auto path = thread_model.path(x_train, y_train, alphas);
 
  if (gamma > 0) {
    path = thread_model.relax(path, x_train, y_train, gamma);
  }
 
  Eigen::ArrayXd scores = Eigen::ArrayXd::Zero(path.size());
 
  for (int j = 0; j < path.size(); ++j) {
    auto eta = path(j).predict(x_test, "linear");
    scores(j) = scorer->eval(eta, y_test, loss);
  }
 
  return scores;
}
 
} // namespace detail
 
template<typename T>
CvResult
crossValidate(Slope model,
              Eigen::EigenBase<T>& x,
              const Eigen::MatrixXd& y_in,
              const CvConfig& config = CvConfig())
{
  CvResult cv_result;
 
  int n = y_in.rows();
 
  auto loss = setupLoss(model.getLossType());
 
  auto y = loss->preprocessResponse(y_in);
  auto scorer = Score::create(config.metric);
 
  auto hyperparams = config.default_hyperparams;
 
  // Override with user-specified parameters
  for (const auto& [key, values] : config.hyperparams) {
    hyperparams[key] = values;
  }
 
  auto grid = detail::createGrid(hyperparams);
 
  // Total number of evaluations (n_repeats * n_folds)
  Folds folds =
    config.predefined_folds.has_value()
      ? Folds(*config.predefined_folds)
      : Folds(n, config.n_folds, config.n_repeats, config.random_seed);
 
  int n_evals = folds.numEvals();
 
  for (const auto& params : grid) {
    GridResult result;
    result.params = params;
 
    double q = params.at("q");
    double gamma = params.at("gamma");
 
    model.setQ(q);
 
    auto initial_path = model.path(x, y);
 
    result.alphas = initial_path.getAlpha();
    int n_alpha = result.alphas.size();
 
    assert((result.alphas > 0).all());
 
    Eigen::MatrixXd scores = Eigen::MatrixXd::Zero(n_evals, n_alpha);
 
#ifdef _OPENMP
    Eigen::setNbThreads(1);
#endif
 
    // Thread-safety for exceptions
    std::vector<std::string> thread_errors(n_evals);
    std::atomic<bool> had_exception{ false };
 
#ifdef _OPENMP
    // omp_set_max_active_levels() is OpenMP 3.0; MSVC's vcomp runtime only
    // implements 2.0, where omp_set_nested() is the way to stop nesting.
#if _OPENMP >= 200805
    omp_set_max_active_levels(1);
#else
    omp_set_nested(0);
#endif
#pragma omp parallel for num_threads(Threads::get())                           \
  shared(scores, thread_errors, had_exception)
#endif
    for (int i = 0; i < n_evals; ++i) {
      try {
        auto [rep, fold] = std::div(i, folds.numFolds());
 
        Slope thread_model = model;
 
        scores.row(i) = detail::fitToFold(x.derived(),
                                          y,
                                          folds,
                                          loss,
                                          scorer,
                                          result.alphas,
                                          thread_model,
                                          fold,
                                          rep,
                                          gamma,
                                          config.copy_x);
 
      } catch (const std::exception& e) {
        thread_errors[i] = e.what();
        had_exception = true;
      } catch (...) {
        thread_errors[i] = "Unknown exception";
        had_exception = true;
      }
    }
 
    if (had_exception) {
      std::string error_message = "Exception(s) during cross-validation:\n";
      for (int i = 0; i < n_evals; ++i) {
        if (!thread_errors[i].empty()) {
          error_message +=
            "Fold " + std::to_string(i) + ": " + thread_errors[i] + "\n";
        }
      }
      throw std::runtime_error(error_message);
    }
 
    result.mean_scores = scores.colwise().mean();
    result.std_errors = stdDevs(scores).array() / std::sqrt(n_evals);
    result.score = std::move(scores);
    cv_result.results.push_back(result);
  }
 
#ifdef _OPENMP
  Eigen::setNbThreads(0);
#endif
 
  detail::findBestParameters(cv_result, scorer);
 
  return cv_result;
}
 
} // namespace slope