cv.h

Go to the documentation of this file.

 
#pragma once
 
#include "folds.h"
#include "score.h"
#include "slope.h"
#include <vector>
 
#ifdef _OPENMP
#include <omp.h>
#endif
 
namespace slope {
 
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;
 
  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;
};
 
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 MatrixType>
CvResult
crossValidate(Slope model,
              MatrixType& 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 = createGrid(hyperparams);
 
  // Total number of evaluations (n_repeats * n_folds)
  Folds folds =
    config.predefined_folds.has_value()
      ? Folds(config.predefined_folds.value())
      : 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);
 
    Eigen::setNbThreads(1);
 
    // Thread-safety for exceptions
    std::vector<std::string> thread_errors(n_evals);
    bool had_exception = false;
 
#ifdef _OPENMP
    omp_set_max_active_levels(1);
#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;
        thread_model.setModifyX(true);
 
        // TODO: Maybe consider not copying at all?
        auto [x_train, y_train, x_test, y_test] = folds.split(x, y, fold, rep);
 
        auto path = thread_model.path(x_train, y_train, result.alphas);
 
        if (gamma > 0) {
          path = thread_model.relax(path, x_train, y_train, gamma);
        }
 
        for (int j = 0; j < n_alpha; ++j) {
          auto eta = path(j).predict(x_test, "linear");
          scores(i, j) = scorer->eval(eta, y_test, loss);
        }
      } catch (const std::exception& e) {
        thread_errors[i] = e.what();
#ifdef _OPENMP
#pragma omp atomic write
#endif
        had_exception = true;
      } catch (...) {
        thread_errors[i] = "Unknown exception";
#ifdef _OPENMP
#pragma omp atomic write
#endif
        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
 
  findBestParameters(cv_result, scorer);
 
  return cv_result;
}
 
} // namespace slope