1// Ceres Solver - A fast non-linear least squares minimizer
2// Copyright 2012 Google Inc. All rights reserved.
3// http://code.google.com/p/ceres-solver/
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
67// In this file, we implement the same evaluation methodology using
68// Ceres. Currently using Levenberg-Marquard with DENSE_QR, we get
69//
70//               Excel  Gnuplot  GaussFit  HBN  MinPack  Ceres
71// Average LRE     2.3      4.3       4.0  6.8      4.4    9.4
72//      Winner       0        0         5   11        2     41
73
74#include <iostream>
75#include <iterator>
76#include <fstream>
77#include "ceres/ceres.h"
78#include "gflags/gflags.h"
79#include "glog/logging.h"
80#include "Eigen/Core"
81
82DEFINE_string(nist_data_dir, "", "Directory containing the NIST non-linear"
83              "regression examples");
84DEFINE_string(minimizer, "trust_region",
85              "Minimizer type to use, choices are: line_search & trust_region");
86DEFINE_string(trust_region_strategy, "levenberg_marquardt",
87              "Options are: levenberg_marquardt, dogleg");
88DEFINE_string(dogleg, "traditional_dogleg",
89              "Options are: traditional_dogleg, subspace_dogleg");
90DEFINE_string(linear_solver, "dense_qr", "Options are: "
91              "sparse_cholesky, dense_qr, dense_normal_cholesky and"
92              "cgnr");
93DEFINE_string(preconditioner, "jacobi", "Options are: "
94              "identity, jacobi");
95DEFINE_string(line_search, "armijo",
96              "Line search algorithm to use, choices are: armijo and wolfe.");
97DEFINE_string(line_search_direction, "lbfgs",
98              "Line search direction algorithm to use, choices: lbfgs, bfgs");
99DEFINE_int32(max_line_search_iterations, 20,
100             "Maximum number of iterations for each line search.");
101DEFINE_int32(max_line_search_restarts, 10,
102             "Maximum number of restarts of line search direction algorithm.");
103DEFINE_string(line_search_interpolation, "cubic",
104              "Degree of polynomial aproximation in line search, "
105              "choices are: bisection, quadratic & cubic.");
106DEFINE_int32(lbfgs_rank, 20,
107             "Rank of L-BFGS inverse Hessian approximation in line search.");
108DEFINE_bool(approximate_eigenvalue_bfgs_scaling, false,
109            "Use approximate eigenvalue scaling in (L)BFGS line search.");
110DEFINE_double(sufficient_decrease, 1.0e-4,
111              "Line search Armijo sufficient (function) decrease factor.");
112DEFINE_double(sufficient_curvature_decrease, 0.9,
113              "Line search Wolfe sufficient curvature decrease factor.");
114DEFINE_int32(num_iterations, 10000, "Number of iterations");
115DEFINE_bool(nonmonotonic_steps, false, "Trust region algorithm can use"
116            " nonmonotic steps");
117DEFINE_double(initial_trust_region_radius, 1e4, "Initial trust region radius");
118
119namespace ceres {
120namespace examples {
121
122using Eigen::Dynamic;
123using Eigen::RowMajor;
124typedef Eigen::Matrix<double, Dynamic, 1> Vector;
125typedef Eigen::Matrix<double, Dynamic, Dynamic, RowMajor> Matrix;
126
127void SplitStringUsingChar(const string& full,
128                          const char delim,
129                          vector<string>* result) {
130  back_insert_iterator< vector<string> > it(*result);
131
132  const char* p = full.data();
133  const char* end = p + full.size();
134  while (p != end) {
135    if (*p == delim) {
136      ++p;
137    } else {
138      const char* start = p;
139      while (++p != end && *p != delim) {
140        // Skip to the next occurence of the delimiter.
141      }
142      *it++ = string(start, p - start);
143    }
144  }
145}
146
147bool GetAndSplitLine(std::ifstream& ifs, std::vector<std::string>* pieces) {
148  pieces->clear();
149  char buf[256];
150  ifs.getline(buf, 256);
151  SplitStringUsingChar(std::string(buf), ' ', pieces);
152  return true;
153}
154
155void SkipLines(std::ifstream& ifs, int num_lines) {
156  char buf[256];
157  for (int i = 0; i < num_lines; ++i) {
158    ifs.getline(buf, 256);
159  }
160}
161
162bool IsSuccessfulTermination(ceres::SolverTerminationType status) {
163  return
164      (status == ceres::FUNCTION_TOLERANCE) ||
165      (status == ceres::GRADIENT_TOLERANCE) ||
166      (status == ceres::PARAMETER_TOLERANCE) ||
167      (status == ceres::USER_SUCCESS);
168}
169
170class NISTProblem {
171 public:
172  explicit NISTProblem(const std::string& filename) {
173    std::ifstream ifs(filename.c_str(), std::ifstream::in);
174
175    std::vector<std::string> pieces;
176    SkipLines(ifs, 24);
177    GetAndSplitLine(ifs, &pieces);
178    const int kNumResponses = std::atoi(pieces[1].c_str());
179
180    GetAndSplitLine(ifs, &pieces);
181    const int kNumPredictors = std::atoi(pieces[0].c_str());
182
183    GetAndSplitLine(ifs, &pieces);
184    const int kNumObservations = std::atoi(pieces[0].c_str());
185
186    SkipLines(ifs, 4);
187    GetAndSplitLine(ifs, &pieces);
188    const int kNumParameters = std::atoi(pieces[0].c_str());
189    SkipLines(ifs, 8);
190
191    // Get the first line of initial and final parameter values to
192    // determine the number of tries.
193    GetAndSplitLine(ifs, &pieces);
194    const int kNumTries = pieces.size() - 4;
195
196    predictor_.resize(kNumObservations, kNumPredictors);
197    response_.resize(kNumObservations, kNumResponses);
198    initial_parameters_.resize(kNumTries, kNumParameters);
199    final_parameters_.resize(1, kNumParameters);
200
201    // Parse the line for parameter b1.
202    int parameter_id = 0;
203    for (int i = 0; i < kNumTries; ++i) {
204      initial_parameters_(i, parameter_id) = std::atof(pieces[i + 2].c_str());
205    }
206    final_parameters_(0, parameter_id) = std::atof(pieces[2 + kNumTries].c_str());
207
208    // Parse the remaining parameter lines.
209    for (int parameter_id = 1; parameter_id < kNumParameters; ++parameter_id) {
210     GetAndSplitLine(ifs, &pieces);
211     // b2, b3, ....
212     for (int i = 0; i < kNumTries; ++i) {
213       initial_parameters_(i, parameter_id) = std::atof(pieces[i + 2].c_str());
214     }
215     final_parameters_(0, parameter_id) = std::atof(pieces[2 + kNumTries].c_str());
216    }
217
218    // Certfied cost
219    SkipLines(ifs, 1);
220    GetAndSplitLine(ifs, &pieces);
221    certified_cost_ = std::atof(pieces[4].c_str()) / 2.0;
222
223    // Read the observations.
224    SkipLines(ifs, 18 - kNumParameters);
225    for (int i = 0; i < kNumObservations; ++i) {
226      GetAndSplitLine(ifs, &pieces);
227      // Response.
228      for (int j = 0; j < kNumResponses; ++j) {
229        response_(i, j) =  std::atof(pieces[j].c_str());
230      }
231
232      // Predictor variables.
233      for (int j = 0; j < kNumPredictors; ++j) {
234        predictor_(i, j) =  std::atof(pieces[j + kNumResponses].c_str());
235      }
236    }
237  }
238
239  Matrix initial_parameters(int start) const { return initial_parameters_.row(start); }
240  Matrix final_parameters() const  { return final_parameters_; }
241  Matrix predictor()        const { return predictor_;         }
242  Matrix response()         const { return response_;          }
243  int predictor_size()      const { return predictor_.cols();  }
244  int num_observations()    const { return predictor_.rows();  }
245  int response_size()       const { return response_.cols();   }
246  int num_parameters()      const { return initial_parameters_.cols(); }
247  int num_starts()          const { return initial_parameters_.rows(); }
248  double certified_cost()   const { return certified_cost_; }
249
250 private:
251  Matrix predictor_;
252  Matrix response_;
253  Matrix initial_parameters_;
254  Matrix final_parameters_;
255  double certified_cost_;
256};
257
258#define NIST_BEGIN(CostFunctionName) \
259  struct CostFunctionName { \
260    CostFunctionName(const double* const x, \
261                     const double* const y) \
262        : x_(*x), y_(*y) {} \
263    double x_; \
264    double y_; \
265    template <typename T> \
266    bool operator()(const T* const b, T* residual) const { \
267    const T y(y_); \
268    const T x(x_); \
269      residual[0] = y - (
270
271#define NIST_END ); return true; }};
272
273// y = b1 * (b2+x)**(-1/b3)  +  e
274NIST_BEGIN(Bennet5)
275  b[0] * pow(b[1] + x, T(-1.0) / b[2])
276NIST_END
277
278// y = b1*(1-exp[-b2*x])  +  e
279NIST_BEGIN(BoxBOD)
280  b[0] * (T(1.0) - exp(-b[1] * x))
281NIST_END
282
283// y = exp[-b1*x]/(b2+b3*x)  +  e
284NIST_BEGIN(Chwirut)
285  exp(-b[0] * x) / (b[1] + b[2] * x)
286NIST_END
287
288// y  = b1*x**b2  +  e
289NIST_BEGIN(DanWood)
290  b[0] * pow(x, b[1])
291NIST_END
292
293// y = b1*exp( -b2*x ) + b3*exp( -(x-b4)**2 / b5**2 )
294//     + b6*exp( -(x-b7)**2 / b8**2 ) + e
295NIST_BEGIN(Gauss)
296  b[0] * exp(-b[1] * x) +
297  b[2] * exp(-pow((x - b[3])/b[4], 2)) +
298  b[5] * exp(-pow((x - b[6])/b[7],2))
299NIST_END
300
301// y = b1*exp(-b2*x) + b3*exp(-b4*x) + b5*exp(-b6*x)  +  e
302NIST_BEGIN(Lanczos)
303  b[0] * exp(-b[1] * x) + b[2] * exp(-b[3] * x) + b[4] * exp(-b[5] * x)
304NIST_END
305
306// y = (b1+b2*x+b3*x**2+b4*x**3) /
307//     (1+b5*x+b6*x**2+b7*x**3)  +  e
308NIST_BEGIN(Hahn1)
309  (b[0] + b[1] * x + b[2] * x * x + b[3] * x * x * x) /
310  (T(1.0) + b[4] * x + b[5] * x * x + b[6] * x * x * x)
311NIST_END
312
313// y = (b1 + b2*x + b3*x**2) /
314//    (1 + b4*x + b5*x**2)  +  e
315NIST_BEGIN(Kirby2)
316  (b[0] + b[1] * x + b[2] * x * x) /
317  (T(1.0) + b[3] * x + b[4] * x * x)
318NIST_END
319
320// y = b1*(x**2+x*b2) / (x**2+x*b3+b4)  +  e
321NIST_BEGIN(MGH09)
322  b[0] * (x * x + x * b[1]) / (x * x + x * b[2] + b[3])
323NIST_END
324
325// y = b1 * exp[b2/(x+b3)]  +  e
326NIST_BEGIN(MGH10)
327  b[0] * exp(b[1] / (x + b[2]))
328NIST_END
329
330// y = b1 + b2*exp[-x*b4] + b3*exp[-x*b5]
331NIST_BEGIN(MGH17)
332  b[0] + b[1] * exp(-x * b[3]) + b[2] * exp(-x * b[4])
333NIST_END
334
335// y = b1*(1-exp[-b2*x])  +  e
336NIST_BEGIN(Misra1a)
337  b[0] * (T(1.0) - exp(-b[1] * x))
338NIST_END
339
340// y = b1 * (1-(1+b2*x/2)**(-2))  +  e
341NIST_BEGIN(Misra1b)
342  b[0] * (T(1.0) - T(1.0)/ ((T(1.0) + b[1] * x / 2.0) * (T(1.0) + b[1] * x / 2.0)))
343NIST_END
344
345// y = b1 * (1-(1+2*b2*x)**(-.5))  +  e
346NIST_BEGIN(Misra1c)
347  b[0] * (T(1.0) - pow(T(1.0) + T(2.0) * b[1] * x, -0.5))
348NIST_END
349
350// y = b1*b2*x*((1+b2*x)**(-1))  +  e
351NIST_BEGIN(Misra1d)
352  b[0] * b[1] * x / (T(1.0) + b[1] * x)
353NIST_END
354
355const double kPi = 3.141592653589793238462643383279;
356// pi = 3.141592653589793238462643383279E0
357// y =  b1 - b2*x - arctan[b3/(x-b4)]/pi  +  e
358NIST_BEGIN(Roszman1)
359  b[0] - b[1] * x - atan2(b[2], (x - b[3]))/T(kPi)
360NIST_END
361
362// y = b1 / (1+exp[b2-b3*x])  +  e
363NIST_BEGIN(Rat42)
364  b[0] / (T(1.0) + exp(b[1] - b[2] * x))
365NIST_END
366
367// y = b1 / ((1+exp[b2-b3*x])**(1/b4))  +  e
368NIST_BEGIN(Rat43)
369  b[0] / pow(T(1.0) + exp(b[1] - b[2] * x), T(1.0) / b[3])
370NIST_END
371
372// y = (b1 + b2*x + b3*x**2 + b4*x**3) /
373//    (1 + b5*x + b6*x**2 + b7*x**3)  +  e
374NIST_BEGIN(Thurber)
375  (b[0] + b[1] * x + b[2] * x * x  + b[3] * x * x * x) /
376  (T(1.0) + b[4] * x + b[5] * x * x + b[6] * x * x * x)
377NIST_END
378
379// y = b1 + b2*cos( 2*pi*x/12 ) + b3*sin( 2*pi*x/12 )
380//        + b5*cos( 2*pi*x/b4 ) + b6*sin( 2*pi*x/b4 )
381//        + b8*cos( 2*pi*x/b7 ) + b9*sin( 2*pi*x/b7 )  + e
382NIST_BEGIN(ENSO)
383  b[0] + b[1] * cos(T(2.0 * kPi) * x / T(12.0)) +
384         b[2] * sin(T(2.0 * kPi) * x / T(12.0)) +
385         b[4] * cos(T(2.0 * kPi) * x / b[3]) +
386         b[5] * sin(T(2.0 * kPi) * x / b[3]) +
387         b[7] * cos(T(2.0 * kPi) * x / b[6]) +
388         b[8] * sin(T(2.0 * kPi) * x / b[6])
389NIST_END
390
391// y = (b1/b2) * exp[-0.5*((x-b3)/b2)**2]  +  e
392NIST_BEGIN(Eckerle4)
393  b[0] / b[1] * exp(T(-0.5) * pow((x - b[2])/b[1], 2))
394NIST_END
395
396struct Nelson {
397 public:
398  Nelson(const double* const x, const double* const y)
399      : x1_(x[0]), x2_(x[1]), y_(y[0]) {}
400
401  template <typename T>
402  bool operator()(const T* const b, T* residual) const {
403    // log[y] = b1 - b2*x1 * exp[-b3*x2]  +  e
404    residual[0] = T(log(y_)) - (b[0] - b[1] * T(x1_) * exp(-b[2] * T(x2_)));
405    return true;
406  }
407
408 private:
409  double x1_;
410  double x2_;
411  double y_;
412};
413
414template <typename Model, int num_residuals, int num_parameters>
415int RegressionDriver(const std::string& filename,
416                     const ceres::Solver::Options& options) {
417  NISTProblem nist_problem(FLAGS_nist_data_dir + filename);
418  CHECK_EQ(num_residuals, nist_problem.response_size());
419  CHECK_EQ(num_parameters, nist_problem.num_parameters());
420
421  Matrix predictor = nist_problem.predictor();
422  Matrix response = nist_problem.response();
423  Matrix final_parameters = nist_problem.final_parameters();
424
425  printf("%s\n", filename.c_str());
426
427  // Each NIST problem comes with multiple starting points, so we
428  // construct the problem from scratch for each case and solve it.
429  int num_success = 0;
430  for (int start = 0; start < nist_problem.num_starts(); ++start) {
431    Matrix initial_parameters = nist_problem.initial_parameters(start);
432
433    ceres::Problem problem;
434    for (int i = 0; i < nist_problem.num_observations(); ++i) {
435      problem.AddResidualBlock(
436          new ceres::AutoDiffCostFunction<Model, num_residuals, num_parameters>(
437              new Model(predictor.data() + nist_problem.predictor_size() * i,
438                        response.data() + nist_problem.response_size() * i)),
439          NULL,
440          initial_parameters.data());
441    }
442
443    ceres::Solver::Summary summary;
444    Solve(options, &problem, &summary);
445
446    // Compute the LRE by comparing each component of the solution
447    // with the ground truth, and taking the minimum.
448    Matrix final_parameters = nist_problem.final_parameters();
449    const double kMaxNumSignificantDigits = 11;
450    double log_relative_error = kMaxNumSignificantDigits + 1;
451    for (int i = 0; i < num_parameters; ++i) {
452      const double tmp_lre =
453          -std::log10(std::fabs(final_parameters(i) - initial_parameters(i)) /
454                      std::fabs(final_parameters(i)));
455      // The maximum LRE is capped at 11 - the precision at which the
456      // ground truth is known.
457      //
458      // The minimum LRE is capped at 0 - no digits match between the
459      // computed solution and the ground truth.
460      log_relative_error =
461          std::min(log_relative_error,
462                   std::max(0.0, std::min(kMaxNumSignificantDigits, tmp_lre)));
463    }
464
465    const int kMinNumMatchingDigits = 4;
466    if (log_relative_error >= kMinNumMatchingDigits) {
467      ++num_success;
468    }
469
470    printf("start: %d status: %s lre: %4.1f initial cost: %e final cost:%e "
471           "certified cost: %e total iterations: %d\n",
472           start + 1,
473           log_relative_error < kMinNumMatchingDigits ? "FAILURE" : "SUCCESS",
474           log_relative_error,
475           summary.initial_cost,
476           summary.final_cost,
477           nist_problem.certified_cost(),
478           (summary.num_successful_steps + summary.num_unsuccessful_steps));
479  }
480  return num_success;
481}
482
483void SetMinimizerOptions(ceres::Solver::Options* options) {
484  CHECK(ceres::StringToMinimizerType(FLAGS_minimizer,
485                                     &options->minimizer_type));
486  CHECK(ceres::StringToLinearSolverType(FLAGS_linear_solver,
487                                        &options->linear_solver_type));
488  CHECK(ceres::StringToPreconditionerType(FLAGS_preconditioner,
489                                          &options->preconditioner_type));
490  CHECK(ceres::StringToTrustRegionStrategyType(
491            FLAGS_trust_region_strategy,
492            &options->trust_region_strategy_type));
493  CHECK(ceres::StringToDoglegType(FLAGS_dogleg, &options->dogleg_type));
494  CHECK(ceres::StringToLineSearchDirectionType(
495      FLAGS_line_search_direction,
496      &options->line_search_direction_type));
497  CHECK(ceres::StringToLineSearchType(FLAGS_line_search,
498                                      &options->line_search_type));
499  CHECK(ceres::StringToLineSearchInterpolationType(
500      FLAGS_line_search_interpolation,
501      &options->line_search_interpolation_type));
502
503  options->max_num_iterations = FLAGS_num_iterations;
504  options->use_nonmonotonic_steps = FLAGS_nonmonotonic_steps;
505  options->initial_trust_region_radius = FLAGS_initial_trust_region_radius;
506  options->max_lbfgs_rank = FLAGS_lbfgs_rank;
507  options->line_search_sufficient_function_decrease = FLAGS_sufficient_decrease;
508  options->line_search_sufficient_curvature_decrease =
509      FLAGS_sufficient_curvature_decrease;
510  options->max_num_line_search_step_size_iterations =
511      FLAGS_max_line_search_iterations;
512  options->max_num_line_search_direction_restarts =
513      FLAGS_max_line_search_restarts;
514  options->use_approximate_eigenvalue_bfgs_scaling =
515      FLAGS_approximate_eigenvalue_bfgs_scaling;
516  options->function_tolerance = 1e-18;
517  options->gradient_tolerance = 1e-18;
518  options->parameter_tolerance = 1e-18;
519}
520
521void SolveNISTProblems() {
522  if (FLAGS_nist_data_dir.empty()) {
523    LOG(FATAL) << "Must specify the directory containing the NIST problems";
524  }
525
526  ceres::Solver::Options options;
527  SetMinimizerOptions(&options);
528
529  std::cout << "Lower Difficulty\n";
530  int easy_success = 0;
531  easy_success += RegressionDriver<Misra1a,  1, 2>("Misra1a.dat",  options);
532  easy_success += RegressionDriver<Chwirut,  1, 3>("Chwirut1.dat", options);
533  easy_success += RegressionDriver<Chwirut,  1, 3>("Chwirut2.dat", options);
534  easy_success += RegressionDriver<Lanczos,  1, 6>("Lanczos3.dat", options);
535  easy_success += RegressionDriver<Gauss,    1, 8>("Gauss1.dat",   options);
536  easy_success += RegressionDriver<Gauss,    1, 8>("Gauss2.dat",   options);
537  easy_success += RegressionDriver<DanWood,  1, 2>("DanWood.dat",  options);
538  easy_success += RegressionDriver<Misra1b,  1, 2>("Misra1b.dat",  options);
539
540  std::cout << "\nMedium Difficulty\n";
541  int medium_success = 0;
542  medium_success += RegressionDriver<Kirby2,   1, 5>("Kirby2.dat",   options);
543  medium_success += RegressionDriver<Hahn1,    1, 7>("Hahn1.dat",    options);
544  medium_success += RegressionDriver<Nelson,   1, 3>("Nelson.dat",   options);
545  medium_success += RegressionDriver<MGH17,    1, 5>("MGH17.dat",    options);
546  medium_success += RegressionDriver<Lanczos,  1, 6>("Lanczos1.dat", options);
547  medium_success += RegressionDriver<Lanczos,  1, 6>("Lanczos2.dat", options);
548  medium_success += RegressionDriver<Gauss,    1, 8>("Gauss3.dat",   options);
549  medium_success += RegressionDriver<Misra1c,  1, 2>("Misra1c.dat",  options);
550  medium_success += RegressionDriver<Misra1d,  1, 2>("Misra1d.dat",  options);
551  medium_success += RegressionDriver<Roszman1, 1, 4>("Roszman1.dat", options);
552  medium_success += RegressionDriver<ENSO,     1, 9>("ENSO.dat",     options);
553
554  std::cout << "\nHigher Difficulty\n";
555  int hard_success = 0;
556  hard_success += RegressionDriver<MGH09,    1, 4>("MGH09.dat",    options);
557  hard_success += RegressionDriver<Thurber,  1, 7>("Thurber.dat",  options);
558  hard_success += RegressionDriver<BoxBOD,   1, 2>("BoxBOD.dat",   options);
559  hard_success += RegressionDriver<Rat42,    1, 3>("Rat42.dat",    options);
560  hard_success += RegressionDriver<MGH10,    1, 3>("MGH10.dat",    options);
561
562  hard_success += RegressionDriver<Eckerle4, 1, 3>("Eckerle4.dat", options);
563  hard_success += RegressionDriver<Rat43,    1, 4>("Rat43.dat",    options);
564  hard_success += RegressionDriver<Bennet5,  1, 3>("Bennett5.dat", options);
565
566  std::cout << "\n";
567  std::cout << "Easy    : " << easy_success << "/16\n";
568  std::cout << "Medium  : " << medium_success << "/22\n";
569  std::cout << "Hard    : " << hard_success << "/16\n";
570  std::cout << "Total   : " << easy_success + medium_success + hard_success << "/54\n";
571}
572
573}  // namespace examples
574}  // namespace ceres
575
576int main(int argc, char** argv) {
577  google::ParseCommandLineFlags(&argc, &argv, true);
578  google::InitGoogleLogging(argv[0]);
579  ceres::examples::SolveNISTProblems();
580  return 0;
581};
582