GoogleGit

blob: f7bd38e1fe341ce6803cc5fbfeb0516a9dd10fa2 [file] [log] [blame]
  1. // Ceres Solver - A fast non-linear least squares minimizer
  2. // Copyright 2015 Google Inc. All rights reserved.
  3. // http://ceres-solver.org/
  4. //
  5. // Redistribution and use in source and binary forms, with or without
  6. // modification, are permitted provided that the following conditions are met:
  7. //
  8. // * Redistributions of source code must retain the above copyright notice,
  9. // this list of conditions and the following disclaimer.
  10. // * Redistributions in binary form must reproduce the above copyright notice,
  11. // this list of conditions and the following disclaimer in the documentation
  12. // and/or other materials provided with the distribution.
  13. // * Neither the name of Google Inc. nor the names of its contributors may be
  14. // used to endorse or promote products derived from this software without
  15. // specific prior written permission.
  16. //
  17. // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  18. // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  19. // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  20. // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  21. // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  22. // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  23. // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  24. // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  25. // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  26. // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  27. // POSSIBILITY OF SUCH DAMAGE.
  28. //
  29. // Author: sameeragarwal@google.com (Sameer Agarwal)
  30. //
  31. // The National Institute of Standards and Technology has released a
  32. // set of problems to test non-linear least squares solvers.
  33. //
  34. // More information about the background on these problems and
  35. // suggested evaluation methodology can be found at:
  36. //
  37. // http://www.itl.nist.gov/div898/strd/nls/nls_info.shtml
  38. //
  39. // The problem data themselves can be found at
  40. //
  41. // http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml
  42. //
  43. // The problems are divided into three levels of difficulty, Easy,
  44. // Medium and Hard. For each problem there are two starting guesses,
  45. // the first one far away from the global minimum and the second
  46. // closer to it.
  47. //
  48. // A problem is considered successfully solved, if every components of
  49. // the solution matches the globally optimal solution in at least 4
  50. // digits or more.
  51. //
  52. // This dataset was used for an evaluation of Non-linear least squares
  53. // solvers:
  54. //
  55. // P. F. Mondragon & B. Borchers, A Comparison of Nonlinear Regression
  56. // Codes, Journal of Modern Applied Statistical Methods, 4(1):343-351,
  57. // 2005.
  58. //
  59. // The results from Mondragon & Borchers can be summarized as
  60. // Excel Gnuplot GaussFit HBN MinPack
  61. // Average LRE 2.3 4.3 4.0 6.8 4.4
  62. // Winner 1 5 12 29 12
  63. //
  64. // Where the row Winner counts, the number of problems for which the
  65. // solver had the highest LRE.
  66. // In this file, we implement the same evaluation methodology using
  67. // Ceres. Currently using Levenberg-Marquardt with DENSE_QR, we get
  68. //
  69. // Excel Gnuplot GaussFit HBN MinPack Ceres
  70. // Average LRE 2.3 4.3 4.0 6.8 4.4 9.4
  71. // Winner 0 0 5 11 2 41
  72. #include <iostream>
  73. #include <iterator>
  74. #include <fstream>
  75. #include "ceres/ceres.h"
  76. #include "gflags/gflags.h"
  77. #include "glog/logging.h"
  78. #include "Eigen/Core"
  79. DEFINE_string(nist_data_dir, "", "Directory containing the NIST non-linear"
  80. "regression examples");
  81. DEFINE_string(minimizer, "trust_region",
  82. "Minimizer type to use, choices are: line_search & trust_region");
  83. DEFINE_string(trust_region_strategy, "levenberg_marquardt",
  84. "Options are: levenberg_marquardt, dogleg");
  85. DEFINE_string(dogleg, "traditional_dogleg",
  86. "Options are: traditional_dogleg, subspace_dogleg");
  87. DEFINE_string(linear_solver, "dense_qr", "Options are: "
  88. "sparse_cholesky, dense_qr, dense_normal_cholesky and"
  89. "cgnr");
  90. DEFINE_string(preconditioner, "jacobi", "Options are: "
  91. "identity, jacobi");
  92. DEFINE_string(line_search, "wolfe",
  93. "Line search algorithm to use, choices are: armijo and wolfe.");
  94. DEFINE_string(line_search_direction, "lbfgs",
  95. "Line search direction algorithm to use, choices: lbfgs, bfgs");
  96. DEFINE_int32(max_line_search_iterations, 20,
  97. "Maximum number of iterations for each line search.");
  98. DEFINE_int32(max_line_search_restarts, 10,
  99. "Maximum number of restarts of line search direction algorithm.");
  100. DEFINE_string(line_search_interpolation, "cubic",
  101. "Degree of polynomial aproximation in line search, "
  102. "choices are: bisection, quadratic & cubic.");
  103. DEFINE_int32(lbfgs_rank, 20,
  104. "Rank of L-BFGS inverse Hessian approximation in line search.");
  105. DEFINE_bool(approximate_eigenvalue_bfgs_scaling, false,
  106. "Use approximate eigenvalue scaling in (L)BFGS line search.");
  107. DEFINE_double(sufficient_decrease, 1.0e-4,
  108. "Line search Armijo sufficient (function) decrease factor.");
  109. DEFINE_double(sufficient_curvature_decrease, 0.9,
  110. "Line search Wolfe sufficient curvature decrease factor.");
  111. DEFINE_int32(num_iterations, 10000, "Number of iterations");
  112. DEFINE_bool(nonmonotonic_steps, false, "Trust region algorithm can use"
  113. " nonmonotic steps");
  114. DEFINE_double(initial_trust_region_radius, 1e4, "Initial trust region radius");
  115. DEFINE_bool(use_numeric_diff, false,
  116. "Use numeric differentiation instead of automatic "
  117. "differentiation.");
  118. DEFINE_string(numeric_diff_method, "ridders", "When using numeric "
  119. "differentiation, selects algorithm. Options are: central, "
  120. "forward, ridders.");
  121. DEFINE_double(ridders_step_size, 1e-9, "Initial step size for Ridders "
  122. "numeric differentiation.");
  123. DEFINE_int32(ridders_extrapolations, 3, "Maximal number of Ridders "
  124. "extrapolations.");
  125. namespace ceres {
  126. namespace examples {
  127. using Eigen::Dynamic;
  128. using Eigen::RowMajor;
  129. typedef Eigen::Matrix<double, Dynamic, 1> Vector;
  130. typedef Eigen::Matrix<double, Dynamic, Dynamic, RowMajor> Matrix;
  131. using std::atof;
  132. using std::atoi;
  133. using std::cout;
  134. using std::ifstream;
  135. using std::string;
  136. using std::vector;
  137. void SplitStringUsingChar(const string& full,
  138. const char delim,
  139. vector<string>* result) {
  140. std::back_insert_iterator< vector<string> > it(*result);
  141. const char* p = full.data();
  142. const char* end = p + full.size();
  143. while (p != end) {
  144. if (*p == delim) {
  145. ++p;
  146. } else {
  147. const char* start = p;
  148. while (++p != end && *p != delim) {
  149. // Skip to the next occurence of the delimiter.
  150. }
  151. *it++ = string(start, p - start);
  152. }
  153. }
  154. }
  155. bool GetAndSplitLine(ifstream& ifs, vector<string>* pieces) {
  156. pieces->clear();
  157. char buf[256];
  158. ifs.getline(buf, 256);
  159. SplitStringUsingChar(string(buf), ' ', pieces);
  160. return true;
  161. }
  162. void SkipLines(ifstream& ifs, int num_lines) {
  163. char buf[256];
  164. for (int i = 0; i < num_lines; ++i) {
  165. ifs.getline(buf, 256);
  166. }
  167. }
  168. class NISTProblem {
  169. public:
  170. explicit NISTProblem(const string& filename) {
  171. ifstream ifs(filename.c_str(), ifstream::in);
  172. vector<string> pieces;
  173. SkipLines(ifs, 24);
  174. GetAndSplitLine(ifs, &pieces);
  175. const int kNumResponses = atoi(pieces[1].c_str());
  176. GetAndSplitLine(ifs, &pieces);
  177. const int kNumPredictors = atoi(pieces[0].c_str());
  178. GetAndSplitLine(ifs, &pieces);
  179. const int kNumObservations = atoi(pieces[0].c_str());
  180. SkipLines(ifs, 4);
  181. GetAndSplitLine(ifs, &pieces);
  182. const int kNumParameters = atoi(pieces[0].c_str());
  183. SkipLines(ifs, 8);
  184. // Get the first line of initial and final parameter values to
  185. // determine the number of tries.
  186. GetAndSplitLine(ifs, &pieces);
  187. const int kNumTries = pieces.size() - 4;
  188. predictor_.resize(kNumObservations, kNumPredictors);
  189. response_.resize(kNumObservations, kNumResponses);
  190. initial_parameters_.resize(kNumTries, kNumParameters);
  191. final_parameters_.resize(1, kNumParameters);
  192. // Parse the line for parameter b1.
  193. int parameter_id = 0;
  194. for (int i = 0; i < kNumTries; ++i) {
  195. initial_parameters_(i, parameter_id) = atof(pieces[i + 2].c_str());
  196. }
  197. final_parameters_(0, parameter_id) = atof(pieces[2 + kNumTries].c_str());
  198. // Parse the remaining parameter lines.
  199. for (int parameter_id = 1; parameter_id < kNumParameters; ++parameter_id) {
  200. GetAndSplitLine(ifs, &pieces);
  201. // b2, b3, ....
  202. for (int i = 0; i < kNumTries; ++i) {
  203. initial_parameters_(i, parameter_id) = atof(pieces[i + 2].c_str());
  204. }
  205. final_parameters_(0, parameter_id) = atof(pieces[2 + kNumTries].c_str());
  206. }
  207. // Certfied cost
  208. SkipLines(ifs, 1);
  209. GetAndSplitLine(ifs, &pieces);
  210. certified_cost_ = atof(pieces[4].c_str()) / 2.0;
  211. // Read the observations.
  212. SkipLines(ifs, 18 - kNumParameters);
  213. for (int i = 0; i < kNumObservations; ++i) {
  214. GetAndSplitLine(ifs, &pieces);
  215. // Response.
  216. for (int j = 0; j < kNumResponses; ++j) {
  217. response_(i, j) = atof(pieces[j].c_str());
  218. }
  219. // Predictor variables.
  220. for (int j = 0; j < kNumPredictors; ++j) {
  221. predictor_(i, j) = atof(pieces[j + kNumResponses].c_str());
  222. }
  223. }
  224. }
  225. Matrix initial_parameters(int start) const { return initial_parameters_.row(start); } // NOLINT
  226. Matrix final_parameters() const { return final_parameters_; }
  227. Matrix predictor() const { return predictor_; }
  228. Matrix response() const { return response_; }
  229. int predictor_size() const { return predictor_.cols(); }
  230. int num_observations() const { return predictor_.rows(); }
  231. int response_size() const { return response_.cols(); }
  232. int num_parameters() const { return initial_parameters_.cols(); }
  233. int num_starts() const { return initial_parameters_.rows(); }
  234. double certified_cost() const { return certified_cost_; }
  235. private:
  236. Matrix predictor_;
  237. Matrix response_;
  238. Matrix initial_parameters_;
  239. Matrix final_parameters_;
  240. double certified_cost_;
  241. };
  242. #define NIST_BEGIN(CostFunctionName) \
  243. struct CostFunctionName { \
  244. CostFunctionName(const double* const x, \
  245. const double* const y) \
  246. : x_(*x), y_(*y) {} \
  247. double x_; \
  248. double y_; \
  249. template <typename T> \
  250. bool operator()(const T* const b, T* residual) const { \
  251. const T y(y_); \
  252. const T x(x_); \
  253. residual[0] = y - (
  254. #define NIST_END ); return true; }};
  255. // y = b1 * (b2+x)**(-1/b3) + e
  256. NIST_BEGIN(Bennet5)
  257. b[0] * pow(b[1] + x, T(-1.0) / b[2])
  258. NIST_END
  259. // y = b1*(1-exp[-b2*x]) + e
  260. NIST_BEGIN(BoxBOD)
  261. b[0] * (T(1.0) - exp(-b[1] * x))
  262. NIST_END
  263. // y = exp[-b1*x]/(b2+b3*x) + e
  264. NIST_BEGIN(Chwirut)
  265. exp(-b[0] * x) / (b[1] + b[2] * x)
  266. NIST_END
  267. // y = b1*x**b2 + e
  268. NIST_BEGIN(DanWood)
  269. b[0] * pow(x, b[1])
  270. NIST_END
  271. // y = b1*exp( -b2*x ) + b3*exp( -(x-b4)**2 / b5**2 )
  272. // + b6*exp( -(x-b7)**2 / b8**2 ) + e
  273. NIST_BEGIN(Gauss)
  274. b[0] * exp(-b[1] * x) +
  275. b[2] * exp(-pow((x - b[3])/b[4], 2)) +
  276. b[5] * exp(-pow((x - b[6])/b[7], 2))
  277. NIST_END
  278. // y = b1*exp(-b2*x) + b3*exp(-b4*x) + b5*exp(-b6*x) + e
  279. NIST_BEGIN(Lanczos)
  280. b[0] * exp(-b[1] * x) + b[2] * exp(-b[3] * x) + b[4] * exp(-b[5] * x)
  281. NIST_END
  282. // y = (b1+b2*x+b3*x**2+b4*x**3) /
  283. // (1+b5*x+b6*x**2+b7*x**3) + e
  284. NIST_BEGIN(Hahn1)
  285. (b[0] + b[1] * x + b[2] * x * x + b[3] * x * x * x) /
  286. (T(1.0) + b[4] * x + b[5] * x * x + b[6] * x * x * x)
  287. NIST_END
  288. // y = (b1 + b2*x + b3*x**2) /
  289. // (1 + b4*x + b5*x**2) + e
  290. NIST_BEGIN(Kirby2)
  291. (b[0] + b[1] * x + b[2] * x * x) /
  292. (T(1.0) + b[3] * x + b[4] * x * x)
  293. NIST_END
  294. // y = b1*(x**2+x*b2) / (x**2+x*b3+b4) + e
  295. NIST_BEGIN(MGH09)
  296. b[0] * (x * x + x * b[1]) / (x * x + x * b[2] + b[3])
  297. NIST_END
  298. // y = b1 * exp[b2/(x+b3)] + e
  299. NIST_BEGIN(MGH10)
  300. b[0] * exp(b[1] / (x + b[2]))
  301. NIST_END
  302. // y = b1 + b2*exp[-x*b4] + b3*exp[-x*b5]
  303. NIST_BEGIN(MGH17)
  304. b[0] + b[1] * exp(-x * b[3]) + b[2] * exp(-x * b[4])
  305. NIST_END
  306. // y = b1*(1-exp[-b2*x]) + e
  307. NIST_BEGIN(Misra1a)
  308. b[0] * (T(1.0) - exp(-b[1] * x))
  309. NIST_END
  310. // y = b1 * (1-(1+b2*x/2)**(-2)) + e
  311. NIST_BEGIN(Misra1b)
  312. b[0] * (T(1.0) - T(1.0)/ ((T(1.0) + b[1] * x / 2.0) * (T(1.0) + b[1] * x / 2.0))) // NOLINT
  313. NIST_END
  314. // y = b1 * (1-(1+2*b2*x)**(-.5)) + e
  315. NIST_BEGIN(Misra1c)
  316. b[0] * (T(1.0) - pow(T(1.0) + T(2.0) * b[1] * x, -0.5))
  317. NIST_END
  318. // y = b1*b2*x*((1+b2*x)**(-1)) + e
  319. NIST_BEGIN(Misra1d)
  320. b[0] * b[1] * x / (T(1.0) + b[1] * x)
  321. NIST_END
  322. const double kPi = 3.141592653589793238462643383279;
  323. // pi = 3.141592653589793238462643383279E0
  324. // y = b1 - b2*x - arctan[b3/(x-b4)]/pi + e
  325. NIST_BEGIN(Roszman1)
  326. b[0] - b[1] * x - atan2(b[2], (x - b[3]))/T(kPi)
  327. NIST_END
  328. // y = b1 / (1+exp[b2-b3*x]) + e
  329. NIST_BEGIN(Rat42)
  330. b[0] / (T(1.0) + exp(b[1] - b[2] * x))
  331. NIST_END
  332. // y = b1 / ((1+exp[b2-b3*x])**(1/b4)) + e
  333. NIST_BEGIN(Rat43)
  334. b[0] / pow(T(1.0) + exp(b[1] - b[2] * x), T(1.0) / b[3])
  335. NIST_END
  336. // y = (b1 + b2*x + b3*x**2 + b4*x**3) /
  337. // (1 + b5*x + b6*x**2 + b7*x**3) + e
  338. NIST_BEGIN(Thurber)
  339. (b[0] + b[1] * x + b[2] * x * x + b[3] * x * x * x) /
  340. (T(1.0) + b[4] * x + b[5] * x * x + b[6] * x * x * x)
  341. NIST_END
  342. // y = b1 + b2*cos( 2*pi*x/12 ) + b3*sin( 2*pi*x/12 )
  343. // + b5*cos( 2*pi*x/b4 ) + b6*sin( 2*pi*x/b4 )
  344. // + b8*cos( 2*pi*x/b7 ) + b9*sin( 2*pi*x/b7 ) + e
  345. NIST_BEGIN(ENSO)
  346. b[0] + b[1] * cos(T(2.0 * kPi) * x / T(12.0)) +
  347. b[2] * sin(T(2.0 * kPi) * x / T(12.0)) +
  348. b[4] * cos(T(2.0 * kPi) * x / b[3]) +
  349. b[5] * sin(T(2.0 * kPi) * x / b[3]) +
  350. b[7] * cos(T(2.0 * kPi) * x / b[6]) +
  351. b[8] * sin(T(2.0 * kPi) * x / b[6])
  352. NIST_END
  353. // y = (b1/b2) * exp[-0.5*((x-b3)/b2)**2] + e
  354. NIST_BEGIN(Eckerle4)
  355. b[0] / b[1] * exp(T(-0.5) * pow((x - b[2])/b[1], 2))
  356. NIST_END
  357. struct Nelson {
  358. public:
  359. Nelson(const double* const x, const double* const y)
  360. : x1_(x[0]), x2_(x[1]), y_(y[0]) {}
  361. template <typename T>
  362. bool operator()(const T* const b, T* residual) const {
  363. // log[y] = b1 - b2*x1 * exp[-b3*x2] + e
  364. residual[0] = T(log(y_)) - (b[0] - b[1] * T(x1_) * exp(-b[2] * T(x2_)));
  365. return true;
  366. }
  367. private:
  368. double x1_;
  369. double x2_;
  370. double y_;
  371. };
  372. static void SetNumericDiffOptions(ceres::NumericDiffOptions* options) {
  373. options->max_num_ridders_extrapolations = FLAGS_ridders_extrapolations;
  374. options->ridders_relative_initial_step_size = FLAGS_ridders_step_size;
  375. }
  376. template <typename Model, int num_residuals, int num_parameters>
  377. int RegressionDriver(const string& filename,
  378. const ceres::Solver::Options& options) {
  379. NISTProblem nist_problem(FLAGS_nist_data_dir + filename);
  380. CHECK_EQ(num_residuals, nist_problem.response_size());
  381. CHECK_EQ(num_parameters, nist_problem.num_parameters());
  382. Matrix predictor = nist_problem.predictor();
  383. Matrix response = nist_problem.response();
  384. Matrix final_parameters = nist_problem.final_parameters();
  385. printf("%s\n", filename.c_str());
  386. // Each NIST problem comes with multiple starting points, so we
  387. // construct the problem from scratch for each case and solve it.
  388. int num_success = 0;
  389. for (int start = 0; start < nist_problem.num_starts(); ++start) {
  390. Matrix initial_parameters = nist_problem.initial_parameters(start);
  391. ceres::Problem problem;
  392. for (int i = 0; i < nist_problem.num_observations(); ++i) {
  393. Model* model = new Model(
  394. predictor.data() + nist_problem.predictor_size() * i,
  395. response.data() + nist_problem.response_size() * i);
  396. ceres::CostFunction* cost_function = NULL;
  397. if (FLAGS_use_numeric_diff) {
  398. ceres::NumericDiffOptions options;
  399. SetNumericDiffOptions(&options);
  400. if (FLAGS_numeric_diff_method == "central") {
  401. cost_function = new NumericDiffCostFunction<Model,
  402. ceres::CENTRAL,
  403. num_residuals,
  404. num_parameters>(
  405. model, ceres::TAKE_OWNERSHIP, num_residuals, options);
  406. } else if (FLAGS_numeric_diff_method == "forward") {
  407. cost_function = new NumericDiffCostFunction<Model,
  408. ceres::FORWARD,
  409. num_residuals,
  410. num_parameters>(
  411. model, ceres::TAKE_OWNERSHIP, num_residuals, options);
  412. } else if (FLAGS_numeric_diff_method == "ridders") {
  413. cost_function = new NumericDiffCostFunction<Model,
  414. ceres::RIDDERS,
  415. num_residuals,
  416. num_parameters>(
  417. model, ceres::TAKE_OWNERSHIP, num_residuals, options);
  418. } else {
  419. LOG(ERROR) << "Invalid numeric diff method specified";
  420. return 0;
  421. }
  422. } else {
  423. cost_function =
  424. new ceres::AutoDiffCostFunction<Model,
  425. num_residuals,
  426. num_parameters>(model);
  427. }
  428. problem.AddResidualBlock(cost_function,
  429. NULL,
  430. initial_parameters.data());
  431. }
  432. ceres::Solver::Summary summary;
  433. Solve(options, &problem, &summary);
  434. // Compute the LRE by comparing each component of the solution
  435. // with the ground truth, and taking the minimum.
  436. Matrix final_parameters = nist_problem.final_parameters();
  437. const double kMaxNumSignificantDigits = 11;
  438. double log_relative_error = kMaxNumSignificantDigits + 1;
  439. for (int i = 0; i < num_parameters; ++i) {
  440. const double tmp_lre =
  441. -std::log10(std::fabs(final_parameters(i) - initial_parameters(i)) /
  442. std::fabs(final_parameters(i)));
  443. // The maximum LRE is capped at 11 - the precision at which the
  444. // ground truth is known.
  445. //
  446. // The minimum LRE is capped at 0 - no digits match between the
  447. // computed solution and the ground truth.
  448. log_relative_error =
  449. std::min(log_relative_error,
  450. std::max(0.0, std::min(kMaxNumSignificantDigits, tmp_lre)));
  451. }
  452. const int kMinNumMatchingDigits = 4;
  453. if (log_relative_error >= kMinNumMatchingDigits) {
  454. ++num_success;
  455. }
  456. printf("start: %d status: %s lre: %4.1f initial cost: %e final cost:%e "
  457. "certified cost: %e total iterations: %d\n",
  458. start + 1,
  459. log_relative_error < kMinNumMatchingDigits ? "FAILURE" : "SUCCESS",
  460. log_relative_error,
  461. summary.initial_cost,
  462. summary.final_cost,
  463. nist_problem.certified_cost(),
  464. (summary.num_successful_steps + summary.num_unsuccessful_steps));
  465. }
  466. return num_success;
  467. }
  468. void SetMinimizerOptions(ceres::Solver::Options* options) {
  469. CHECK(ceres::StringToMinimizerType(FLAGS_minimizer,
  470. &options->minimizer_type));
  471. CHECK(ceres::StringToLinearSolverType(FLAGS_linear_solver,
  472. &options->linear_solver_type));
  473. CHECK(ceres::StringToPreconditionerType(FLAGS_preconditioner,
  474. &options->preconditioner_type));
  475. CHECK(ceres::StringToTrustRegionStrategyType(
  476. FLAGS_trust_region_strategy,
  477. &options->trust_region_strategy_type));
  478. CHECK(ceres::StringToDoglegType(FLAGS_dogleg, &options->dogleg_type));
  479. CHECK(ceres::StringToLineSearchDirectionType(
  480. FLAGS_line_search_direction,
  481. &options->line_search_direction_type));
  482. CHECK(ceres::StringToLineSearchType(FLAGS_line_search,
  483. &options->line_search_type));
  484. CHECK(ceres::StringToLineSearchInterpolationType(
  485. FLAGS_line_search_interpolation,
  486. &options->line_search_interpolation_type));
  487. options->max_num_iterations = FLAGS_num_iterations;
  488. options->use_nonmonotonic_steps = FLAGS_nonmonotonic_steps;
  489. options->initial_trust_region_radius = FLAGS_initial_trust_region_radius;
  490. options->max_lbfgs_rank = FLAGS_lbfgs_rank;
  491. options->line_search_sufficient_function_decrease = FLAGS_sufficient_decrease;
  492. options->line_search_sufficient_curvature_decrease =
  493. FLAGS_sufficient_curvature_decrease;
  494. options->max_num_line_search_step_size_iterations =
  495. FLAGS_max_line_search_iterations;
  496. options->max_num_line_search_direction_restarts =
  497. FLAGS_max_line_search_restarts;
  498. options->use_approximate_eigenvalue_bfgs_scaling =
  499. FLAGS_approximate_eigenvalue_bfgs_scaling;
  500. options->function_tolerance = 1e-18;
  501. options->gradient_tolerance = 1e-18;
  502. options->parameter_tolerance = 1e-18;
  503. }
  504. void SolveNISTProblems() {
  505. if (FLAGS_nist_data_dir.empty()) {
  506. LOG(FATAL) << "Must specify the directory containing the NIST problems";
  507. }
  508. ceres::Solver::Options options;
  509. SetMinimizerOptions(&options);
  510. cout << "Lower Difficulty\n";
  511. int easy_success = 0;
  512. easy_success += RegressionDriver<Misra1a, 1, 2>("Misra1a.dat", options);
  513. easy_success += RegressionDriver<Chwirut, 1, 3>("Chwirut1.dat", options);
  514. easy_success += RegressionDriver<Chwirut, 1, 3>("Chwirut2.dat", options);
  515. easy_success += RegressionDriver<Lanczos, 1, 6>("Lanczos3.dat", options);
  516. easy_success += RegressionDriver<Gauss, 1, 8>("Gauss1.dat", options);
  517. easy_success += RegressionDriver<Gauss, 1, 8>("Gauss2.dat", options);
  518. easy_success += RegressionDriver<DanWood, 1, 2>("DanWood.dat", options);
  519. easy_success += RegressionDriver<Misra1b, 1, 2>("Misra1b.dat", options);
  520. cout << "\nMedium Difficulty\n";
  521. int medium_success = 0;
  522. medium_success += RegressionDriver<Kirby2, 1, 5>("Kirby2.dat", options);
  523. medium_success += RegressionDriver<Hahn1, 1, 7>("Hahn1.dat", options);
  524. medium_success += RegressionDriver<Nelson, 1, 3>("Nelson.dat", options);
  525. medium_success += RegressionDriver<MGH17, 1, 5>("MGH17.dat", options);
  526. medium_success += RegressionDriver<Lanczos, 1, 6>("Lanczos1.dat", options);
  527. medium_success += RegressionDriver<Lanczos, 1, 6>("Lanczos2.dat", options);
  528. medium_success += RegressionDriver<Gauss, 1, 8>("Gauss3.dat", options);
  529. medium_success += RegressionDriver<Misra1c, 1, 2>("Misra1c.dat", options);
  530. medium_success += RegressionDriver<Misra1d, 1, 2>("Misra1d.dat", options);
  531. medium_success += RegressionDriver<Roszman1, 1, 4>("Roszman1.dat", options);
  532. medium_success += RegressionDriver<ENSO, 1, 9>("ENSO.dat", options);
  533. cout << "\nHigher Difficulty\n";
  534. int hard_success = 0;
  535. hard_success += RegressionDriver<MGH09, 1, 4>("MGH09.dat", options);
  536. hard_success += RegressionDriver<Thurber, 1, 7>("Thurber.dat", options);
  537. hard_success += RegressionDriver<BoxBOD, 1, 2>("BoxBOD.dat", options);
  538. hard_success += RegressionDriver<Rat42, 1, 3>("Rat42.dat", options);
  539. hard_success += RegressionDriver<MGH10, 1, 3>("MGH10.dat", options);
  540. hard_success += RegressionDriver<Eckerle4, 1, 3>("Eckerle4.dat", options);
  541. hard_success += RegressionDriver<Rat43, 1, 4>("Rat43.dat", options);
  542. hard_success += RegressionDriver<Bennet5, 1, 3>("Bennett5.dat", options);
  543. cout << "\n";
  544. cout << "Easy : " << easy_success << "/16\n";
  545. cout << "Medium : " << medium_success << "/22\n";
  546. cout << "Hard : " << hard_success << "/16\n";
  547. cout << "Total : "
  548. << easy_success + medium_success + hard_success << "/54\n";
  549. }
  550. } // namespace examples
  551. } // namespace ceres
  552. int main(int argc, char** argv) {
  553. CERES_GFLAGS_NAMESPACE::ParseCommandLineFlags(&argc, &argv, true);
  554. google::InitGoogleLogging(argv[0]);
  555. ceres::examples::SolveNISTProblems();
  556. return 0;
  557. }