examples/nist.cc

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