Refactor nist.cc to be compatible with TinySolver Change-Id: Iec0455ff9fe327fe75dc63f5b80c2ecca2c48e55
diff --git a/examples/nist.cc b/examples/nist.cc index 55a3c7c..754467c 100644 --- a/examples/nist.cc +++ b/examples/nist.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. +// Copyright 2017 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -71,14 +71,18 @@ // Average LRE 2.3 4.3 4.0 6.8 4.4 9.4 // Winner 0 0 5 11 2 41 +#include <Eigen/Core> +#include <fstream> #include <iostream> #include <iterator> -#include <fstream> + #include "ceres/ceres.h" +#include "ceres/tiny_solver.h" +#include "ceres/tiny_solver_cost_function_adapter.h" #include "gflags/gflags.h" #include "glog/logging.h" -#include "Eigen/Core" +DEFINE_bool(use_tiny_solver, false, "Use TinySolver instead of Ceres::Solver"); DEFINE_string(nist_data_dir, "", "Directory containing the NIST non-linear" "regression examples"); DEFINE_string(minimizer, "trust_region", @@ -265,20 +269,22 @@ double certified_cost_; }; -#define NIST_BEGIN(CostFunctionName) \ - struct CostFunctionName { \ - CostFunctionName(const double* const x, \ - const double* const y) \ - : x_(*x), y_(*y) {} \ - double x_; \ - double y_; \ - template <typename T> \ - bool operator()(const T* const b, T* residual) const { \ - const T y(y_); \ - const T x(x_); \ - residual[0] = y - ( +#define NIST_BEGIN(CostFunctionName) \ + struct CostFunctionName { \ + CostFunctionName(const double* const x, \ + const double* const y, \ + const int n) \ + : x_(x), y_(y), n_(n) {} \ + const double* x_; \ + const double* y_; \ + const int n_; \ + template <typename T> \ + bool operator()(const T* const b, T* residual) const { \ + for (int i = 0; i < n_; ++i) { \ + const T x(x_[i]); \ + residual[i] = y_[i] - ( -#define NIST_END ); return true; }}; +#define NIST_END ); } return true; }}; // y = b1 * (b2+x)**(-1/b3) + e NIST_BEGIN(Bennet5) @@ -405,20 +411,22 @@ struct Nelson { public: - Nelson(const double* const x, const double* const y) - : x1_(x[0]), x2_(x[1]), y_(y[0]) {} + Nelson(const double* const x, const double* const y, const int n) + : x_(x), y_(y), n_(n) {} template <typename T> bool operator()(const T* const b, T* residual) const { // log[y] = b1 - b2*x1 * exp[-b3*x2] + e - residual[0] = log(y_) - (b[0] - b[1] * x1_ * exp(-b[2] * x2_)); + for (int i = 0; i < n_; ++i) { + residual[i] = log(y_[i]) - (b[0] - b[1] * x_[2 * i] * exp(-b[2] * x_[2 * i + 1])); + } return true; } private: - double x1_; - double x2_; - double y_; + const double* x_; + const double* y_; + const int n_; }; static void SetNumericDiffOptions(ceres::NumericDiffOptions* options) { @@ -426,138 +434,18 @@ options->ridders_relative_initial_step_size = FLAGS_ridders_step_size; } -string JoinPath(const string& dirname, const string& basename) { -#ifdef _WIN32 - static const char separator = '\\'; -#else - static const char separator = '/'; -#endif // _WIN32 - - if ((!basename.empty() && basename[0] == separator) || dirname.empty()) { - return basename; - } else if (dirname[dirname.size() - 1] == separator) { - return dirname + basename; - } else { - return dirname + string(&separator, 1) + basename; - } -} - -template <typename Model, int num_residuals, int num_parameters> -int RegressionDriver(const string& filename, - const ceres::Solver::Options& options) { - NISTProblem nist_problem(JoinPath(FLAGS_nist_data_dir, filename)); - CHECK_EQ(num_residuals, nist_problem.response_size()); - CHECK_EQ(num_parameters, nist_problem.num_parameters()); - - Matrix predictor = nist_problem.predictor(); - Matrix response = nist_problem.response(); - Matrix final_parameters = nist_problem.final_parameters(); - - printf("%s\n", filename.c_str()); - - // Each NIST problem comes with multiple starting points, so we - // construct the problem from scratch for each case and solve it. - int num_success = 0; - for (int start = 0; start < nist_problem.num_starts(); ++start) { - Matrix initial_parameters = nist_problem.initial_parameters(start); - - ceres::Problem problem; - for (int i = 0; i < nist_problem.num_observations(); ++i) { - Model* model = new Model( - predictor.data() + nist_problem.predictor_size() * i, - response.data() + nist_problem.response_size() * i); - ceres::CostFunction* cost_function = NULL; - if (FLAGS_use_numeric_diff) { - ceres::NumericDiffOptions options; - SetNumericDiffOptions(&options); - if (FLAGS_numeric_diff_method == "central") { - cost_function = new NumericDiffCostFunction<Model, - ceres::CENTRAL, - num_residuals, - num_parameters>( - model, ceres::TAKE_OWNERSHIP, num_residuals, options); - } else if (FLAGS_numeric_diff_method == "forward") { - cost_function = new NumericDiffCostFunction<Model, - ceres::FORWARD, - num_residuals, - num_parameters>( - model, ceres::TAKE_OWNERSHIP, num_residuals, options); - } else if (FLAGS_numeric_diff_method == "ridders") { - cost_function = new NumericDiffCostFunction<Model, - ceres::RIDDERS, - num_residuals, - num_parameters>( - model, ceres::TAKE_OWNERSHIP, num_residuals, options); - } else { - LOG(ERROR) << "Invalid numeric diff method specified"; - return 0; - } - } else { - cost_function = - new ceres::AutoDiffCostFunction<Model, - num_residuals, - num_parameters>(model); - } - - problem.AddResidualBlock(cost_function, - NULL, - initial_parameters.data()); - } - - ceres::Solver::Summary summary; - Solve(options, &problem, &summary); - - // Compute the LRE by comparing each component of the solution - // with the ground truth, and taking the minimum. - Matrix final_parameters = nist_problem.final_parameters(); - const double kMaxNumSignificantDigits = 11; - double log_relative_error = kMaxNumSignificantDigits + 1; - for (int i = 0; i < num_parameters; ++i) { - const double tmp_lre = - -std::log10(std::fabs(final_parameters(i) - initial_parameters(i)) / - std::fabs(final_parameters(i))); - // The maximum LRE is capped at 11 - the precision at which the - // ground truth is known. - // - // The minimum LRE is capped at 0 - no digits match between the - // computed solution and the ground truth. - log_relative_error = - std::min(log_relative_error, - std::max(0.0, std::min(kMaxNumSignificantDigits, tmp_lre))); - } - - const int kMinNumMatchingDigits = 4; - if (log_relative_error > kMinNumMatchingDigits) { - ++num_success; - } - - printf("start: %d status: %s lre: %4.1f initial cost: %e final cost:%e " - "certified cost: %e total iterations: %d\n", - start + 1, - log_relative_error < kMinNumMatchingDigits ? "FAILURE" : "SUCCESS", - log_relative_error, - summary.initial_cost, - summary.final_cost, - nist_problem.certified_cost(), - (summary.num_successful_steps + summary.num_unsuccessful_steps)); - } - return num_success; -} - void SetMinimizerOptions(ceres::Solver::Options* options) { - CHECK(ceres::StringToMinimizerType(FLAGS_minimizer, - &options->minimizer_type)); + CHECK( + ceres::StringToMinimizerType(FLAGS_minimizer, &options->minimizer_type)); CHECK(ceres::StringToLinearSolverType(FLAGS_linear_solver, &options->linear_solver_type)); CHECK(ceres::StringToPreconditionerType(FLAGS_preconditioner, &options->preconditioner_type)); CHECK(ceres::StringToTrustRegionStrategyType( - FLAGS_trust_region_strategy, - &options->trust_region_strategy_type)); + FLAGS_trust_region_strategy, &options->trust_region_strategy_type)); CHECK(ceres::StringToDoglegType(FLAGS_dogleg, &options->dogleg_type)); CHECK(ceres::StringToLineSearchDirectionType( - FLAGS_line_search_direction, - &options->line_search_direction_type)); + FLAGS_line_search_direction, &options->line_search_direction_type)); CHECK(ceres::StringToLineSearchType(FLAGS_line_search, &options->line_search_type)); CHECK(ceres::StringToLineSearchInterpolationType( @@ -582,57 +470,213 @@ options->parameter_tolerance = std::numeric_limits<double>::epsilon(); } +string JoinPath(const string& dirname, const string& basename) { +#ifdef _WIN32 + static const char separator = '\\'; +#else + static const char separator = '/'; +#endif // _WIN32 + + if ((!basename.empty() && basename[0] == separator) || dirname.empty()) { + return basename; + } else if (dirname[dirname.size() - 1] == separator) { + return dirname + basename; + } else { + return dirname + string(&separator, 1) + basename; + } +} + +template <typename Model, int num_parameters> +CostFunction* CreateCostFunction(const Matrix& predictor, + const Matrix& response, + const int num_observations) { + Model* model = + new Model(predictor.data(), response.data(), num_observations); + ceres::CostFunction* cost_function = NULL; + if (FLAGS_use_numeric_diff) { + ceres::NumericDiffOptions options; + SetNumericDiffOptions(&options); + if (FLAGS_numeric_diff_method == "central") { + cost_function = new NumericDiffCostFunction<Model, + ceres::CENTRAL, + ceres::DYNAMIC, + num_parameters>( + model, + ceres::TAKE_OWNERSHIP, + num_observations, + options); + } else if (FLAGS_numeric_diff_method == "forward") { + cost_function = new NumericDiffCostFunction<Model, + ceres::FORWARD, + ceres::DYNAMIC, + num_parameters>( + model, + ceres::TAKE_OWNERSHIP, + num_observations, + options); + } else if (FLAGS_numeric_diff_method == "ridders") { + cost_function = new NumericDiffCostFunction<Model, + ceres::RIDDERS, + ceres::DYNAMIC, + num_parameters>( + model, + ceres::TAKE_OWNERSHIP, + num_observations, + options); + } else { + LOG(ERROR) << "Invalid numeric diff method specified"; + return 0; + } + } else { + cost_function = + new ceres::AutoDiffCostFunction<Model, ceres::DYNAMIC, num_parameters>( + model, num_observations); + } + return cost_function; +} + +double ComputeLRE(const Matrix& expected, const Matrix& actual) { + // Compute the LRE by comparing each component of the solution + // with the ground truth, and taking the minimum. + const double kMaxNumSignificantDigits = 11; + double log_relative_error = kMaxNumSignificantDigits + 1; + for (int i = 0; i < expected.cols(); ++i) { + const double tmp_lre = -std::log10(std::fabs(expected(i) - actual(i)) / + std::fabs(expected(i))); + // The maximum LRE is capped at 11 - the precision at which the + // ground truth is known. + // + // The minimum LRE is capped at 0 - no digits match between the + // computed solution and the ground truth. + log_relative_error = + std::min(log_relative_error, + std::max(0.0, std::min(kMaxNumSignificantDigits, tmp_lre))); + } + return log_relative_error; +} + +template <typename Model, int num_parameters> +int RegressionDriver(const string& filename) { + NISTProblem nist_problem(JoinPath(FLAGS_nist_data_dir, filename)); + CHECK_EQ(num_parameters, nist_problem.num_parameters()); + + Matrix predictor = nist_problem.predictor(); + Matrix response = nist_problem.response(); + Matrix final_parameters = nist_problem.final_parameters(); + + printf("%s\n", filename.c_str()); + + // Each NIST problem comes with multiple starting points, so we + // construct the problem from scratch for each case and solve it. + int num_success = 0; + for (int start = 0; start < nist_problem.num_starts(); ++start) { + Matrix initial_parameters = nist_problem.initial_parameters(start); + ceres::CostFunction* cost_function = CreateCostFunction<Model, num_parameters>( + predictor, response, nist_problem.num_observations()); + + double initial_cost; + double final_cost; + + if (!FLAGS_use_tiny_solver) { + ceres::Problem problem; + problem.AddResidualBlock(cost_function, NULL, initial_parameters.data()); + ceres::Solver::Summary summary; + ceres::Solver::Options options; + SetMinimizerOptions(&options); + Solve(options, &problem, &summary); + initial_cost = summary.initial_cost; + final_cost = summary.final_cost; + } else { + ceres::TinySolverCostFunctionAdapter<Eigen::Dynamic, num_parameters> cfa( + *cost_function); + typedef + ceres::TinySolver< + ceres::TinySolverCostFunctionAdapter<Eigen::Dynamic, num_parameters>, + Eigen::LDLT<Eigen::Matrix<double, num_parameters, num_parameters>>> Solver; + Solver solver; + solver.options.max_iterations = FLAGS_num_iterations; + solver.options.error_threshold = std::numeric_limits<double>::epsilon(); + solver.options.gradient_threshold = std::numeric_limits<double>::epsilon(); + solver.options.relative_step_threshold = std::numeric_limits<double>::epsilon(); + + Eigen::Matrix<double, num_parameters,1> x; + x = initial_parameters.transpose(); + typename Solver::Summary summary = solver.Solve(cfa, &x); + initial_parameters = x; + initial_cost = summary.initial_cost; + final_cost = summary.final_cost; + delete cost_function; + } + + const double log_relative_error = ComputeLRE(nist_problem.final_parameters(), + initial_parameters); + const int kMinNumMatchingDigits = 4; + if (log_relative_error > kMinNumMatchingDigits) { + ++num_success; + } + + printf( + "start: %d status: %s lre: %4.1f initial cost: %e final cost:%e " + "certified cost: %e\n", + start + 1, + log_relative_error < kMinNumMatchingDigits ? "FAILURE" : "SUCCESS", + log_relative_error, + initial_cost, + final_cost, + nist_problem.certified_cost()); + } + return num_success; +} + + void SolveNISTProblems() { if (FLAGS_nist_data_dir.empty()) { LOG(FATAL) << "Must specify the directory containing the NIST problems"; } - ceres::Solver::Options options; - SetMinimizerOptions(&options); - cout << "Lower Difficulty\n"; int easy_success = 0; - easy_success += RegressionDriver<Misra1a, 1, 2>("Misra1a.dat", options); - easy_success += RegressionDriver<Chwirut, 1, 3>("Chwirut1.dat", options); - easy_success += RegressionDriver<Chwirut, 1, 3>("Chwirut2.dat", options); - easy_success += RegressionDriver<Lanczos, 1, 6>("Lanczos3.dat", options); - easy_success += RegressionDriver<Gauss, 1, 8>("Gauss1.dat", options); - easy_success += RegressionDriver<Gauss, 1, 8>("Gauss2.dat", options); - easy_success += RegressionDriver<DanWood, 1, 2>("DanWood.dat", options); - easy_success += RegressionDriver<Misra1b, 1, 2>("Misra1b.dat", options); + easy_success += RegressionDriver<Misra1a, 2>("Misra1a.dat"); + easy_success += RegressionDriver<Chwirut, 3>("Chwirut1.dat"); + easy_success += RegressionDriver<Chwirut, 3>("Chwirut2.dat"); + easy_success += RegressionDriver<Lanczos, 6>("Lanczos3.dat"); + easy_success += RegressionDriver<Gauss, 8>("Gauss1.dat"); + easy_success += RegressionDriver<Gauss, 8>("Gauss2.dat"); + easy_success += RegressionDriver<DanWood, 2>("DanWood.dat"); + easy_success += RegressionDriver<Misra1b, 2>("Misra1b.dat"); cout << "\nMedium Difficulty\n"; int medium_success = 0; - medium_success += RegressionDriver<Kirby2, 1, 5>("Kirby2.dat", options); - medium_success += RegressionDriver<Hahn1, 1, 7>("Hahn1.dat", options); - medium_success += RegressionDriver<Nelson, 1, 3>("Nelson.dat", options); - medium_success += RegressionDriver<MGH17, 1, 5>("MGH17.dat", options); - medium_success += RegressionDriver<Lanczos, 1, 6>("Lanczos1.dat", options); - medium_success += RegressionDriver<Lanczos, 1, 6>("Lanczos2.dat", options); - medium_success += RegressionDriver<Gauss, 1, 8>("Gauss3.dat", options); - medium_success += RegressionDriver<Misra1c, 1, 2>("Misra1c.dat", options); - medium_success += RegressionDriver<Misra1d, 1, 2>("Misra1d.dat", options); - medium_success += RegressionDriver<Roszman1, 1, 4>("Roszman1.dat", options); - medium_success += RegressionDriver<ENSO, 1, 9>("ENSO.dat", options); + medium_success += RegressionDriver<Kirby2, 5>("Kirby2.dat"); + medium_success += RegressionDriver<Hahn1, 7>("Hahn1.dat"); + medium_success += RegressionDriver<Nelson, 3>("Nelson.dat"); + medium_success += RegressionDriver<MGH17, 5>("MGH17.dat"); + medium_success += RegressionDriver<Lanczos, 6>("Lanczos1.dat"); + medium_success += RegressionDriver<Lanczos, 6>("Lanczos2.dat"); + medium_success += RegressionDriver<Gauss, 8>("Gauss3.dat"); + medium_success += RegressionDriver<Misra1c, 2>("Misra1c.dat"); + medium_success += RegressionDriver<Misra1d, 2>("Misra1d.dat"); + medium_success += RegressionDriver<Roszman1, 4>("Roszman1.dat"); + medium_success += RegressionDriver<ENSO, 9>("ENSO.dat"); cout << "\nHigher Difficulty\n"; int hard_success = 0; - hard_success += RegressionDriver<MGH09, 1, 4>("MGH09.dat", options); - hard_success += RegressionDriver<Thurber, 1, 7>("Thurber.dat", options); - hard_success += RegressionDriver<BoxBOD, 1, 2>("BoxBOD.dat", options); - hard_success += RegressionDriver<Rat42, 1, 3>("Rat42.dat", options); - hard_success += RegressionDriver<MGH10, 1, 3>("MGH10.dat", options); + hard_success += RegressionDriver<MGH09, 4>("MGH09.dat"); + hard_success += RegressionDriver<Thurber, 7>("Thurber.dat"); + hard_success += RegressionDriver<BoxBOD, 2>("BoxBOD.dat"); + hard_success += RegressionDriver<Rat42, 3>("Rat42.dat"); + hard_success += RegressionDriver<MGH10, 3>("MGH10.dat"); - hard_success += RegressionDriver<Eckerle4, 1, 3>("Eckerle4.dat", options); - hard_success += RegressionDriver<Rat43, 1, 4>("Rat43.dat", options); - hard_success += RegressionDriver<Bennet5, 1, 3>("Bennett5.dat", options); + hard_success += RegressionDriver<Eckerle4, 3>("Eckerle4.dat"); + hard_success += RegressionDriver<Rat43, 4>("Rat43.dat"); + hard_success += RegressionDriver<Bennet5, 3>("Bennett5.dat"); cout << "\n"; cout << "Easy : " << easy_success << "/16\n"; cout << "Medium : " << medium_success << "/22\n"; cout << "Hard : " << hard_success << "/16\n"; - cout << "Total : " - << easy_success + medium_success + hard_success << "/54\n"; + cout << "Total : " << easy_success + medium_success + hard_success + << "/54\n"; } } // namespace examples