cv.h

Go to the documentation of this file.

 
#pragma once
 
#include "eigen_compat.h"
#include "folds.h"
#include "observation_weights.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
scoreFoldPath(const SlopePath& path,
              Eigen::EigenBase<T>& x_test,
              const Eigen::MatrixXd& y_test,
              const std::unique_ptr<Loss>& loss,
              const std::unique_ptr<Score>& scorer,
              const Eigen::VectorXd& observation_weights,
              const bool score_class_labels)
{
  Eigen::MatrixXd targets = y_test;
  if (score_class_labels) {
    targets = loss->predict(loss->link(y_test));
  }
  Eigen::ArrayXd scores(path.size());
  for (int j = 0; j < path.size(); ++j) {
    const auto eta = path(j).predict(x_test, "linear");
    scores(j) = scorer->eval(eta, targets, loss, observation_weights);
  }
  return scores;
}
 
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,
          const bool copy_x,
          const Eigen::VectorXd& observation_weights,
          const bool score_class_labels)
{
  const auto train_idx = folds.getTrainingIndices(fold, rep);
  const auto test_idx = folds.getTestIndices(fold, rep);
  const Eigen::VectorXd train_weights = observation_weights(train_idx);
  const Eigen::VectorXd test_weights = observation_weights(test_idx);
 
  if (copy_x) {
    // Relaxation must reuse the original predictors and their normalization.
    thread_model.setModifyX(gamma == 0);
    auto [x_train, y_train, x_test, y_test] = folds.split(x, y, fold, rep);
    auto path =
      SlopeAccess::path(thread_model, x_train, y_train, alphas, train_weights);
    if (gamma > 0) {
      path = SlopeAccess::relax(
        thread_model, path, x_train, y_train, gamma, train_weights);
    }
    return scoreFoldPath(
      path, x_test, y_test, loss, scorer, test_weights, score_class_labels);
  }
 
  thread_model.setModifyX(false);
  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 =
    SlopeAccess::path(thread_model, x_train, y_train, alphas, train_weights);
  if (gamma > 0) {
    path = SlopeAccess::relax(
      thread_model, path, x_train, y_train, gamma, train_weights);
  }
  return scoreFoldPath(
    path, x_test, y_test, loss, scorer, test_weights, score_class_labels);
}
 
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,
          const bool copy_x,
          const Eigen::VectorXd& observation_weights,
          const bool score_class_labels)
{
  thread_model.setModifyX(gamma == 0);
  const auto train_idx = folds.getTrainingIndices(fold, rep);
  const auto test_idx = folds.getTestIndices(fold, rep);
  const Eigen::VectorXd train_weights = observation_weights(train_idx);
  const Eigen::VectorXd test_weights = observation_weights(test_idx);
  auto [x_train, y_train, x_test, y_test] = folds.split(x, y, fold, rep);
  auto path =
    SlopeAccess::path(thread_model, x_train, y_train, alphas, train_weights);
  if (gamma > 0) {
    path = SlopeAccess::relax(
      thread_model, path, x_train, y_train, gamma, train_weights);
  }
  return scoreFoldPath(
    path, x_test, y_test, loss, scorer, test_weights, score_class_labels);
}
 
} // namespace detail

template<typename T>
CvResult
crossValidate(Slope model,
              Eigen::EigenBase<T>& x,
              const Eigen::MatrixXd& y_in,
              const CvConfig& config = CvConfig(),
              const Eigen::VectorXd& observation_weights = Eigen::VectorXd(0))
{
  CvResult cv_result;
 
  int n = y_in.rows();
  if (x.rows() != n) {
    throw std::invalid_argument("x and y must have the same number of rows");
  }
  if (!isFinite(x.derived())) {
    throw std::invalid_argument("x must not contain NA, NaN, or Inf values");
  }
  if (!y_in.allFinite()) {
    throw std::invalid_argument("y must not contain NA, NaN, or Inf values");
  }
  const Eigen::VectorXd weights =
    detail::normalizeObservationWeights(observation_weights, n);
 
  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 (int rep = 0; rep < folds.numRepetitions(); ++rep) {
    for (int fold = 0; fold < folds.numFolds(); ++fold) {
      const auto train_indices = folds.getTrainingIndices(fold, rep);
      const auto test_indices = folds.getTestIndices(fold, rep);
      const Eigen::VectorXd train_weights = weights(train_indices);
      const Eigen::VectorXd test_weights = weights(test_indices);
      const std::string fold_context = "fold " + std::to_string(fold + 1) +
                                       " in repeat " + std::to_string(rep + 1);
      if (train_weights.sum() <= 0) {
        throw std::invalid_argument("Training " + fold_context +
                                    " has zero total observation weight");
      }
      if (test_weights.sum() <= 0) {
        throw std::invalid_argument("Validation " + fold_context +
                                    " has zero total observation weight");
      }
    }
  }
  const bool score_class_labels =
    config.metric == "accuracy" || config.metric == "misclass";
 
  // Full-data alpha generation must leave predictors available for each fold.
  model.setModifyX(false);
  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 =
      detail::SlopeAccess::path(model, x, y, Eigen::ArrayXd(0), weights);
 
    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;
        // The full-data path has already converted its alphas to absolute
        // values.
        thread_model.setAlphaScale("absolute");
 
        scores.row(i) = detail::fitToFold(x.derived(),
                                          y,
                                          folds,
                                          loss,
                                          scorer,
                                          result.alphas,
                                          thread_model,
                                          fold,
                                          rep,
                                          gamma,
                                          config.copy_x,
                                          weights,
                                          score_class_labels);
 
      } 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