Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 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 | // |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 31 | // The National Institute of Standards and Technology has released a |
| 32 | // set of problems to test non-linear least squares solvers. |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 33 | // |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 34 | // More information about the background on these problems and |
| 35 | // suggested evaluation methodology can be found at: |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 36 | // |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 37 | // http://www.itl.nist.gov/div898/strd/nls/nls_info.shtml |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 38 | // |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 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 |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 73 | |
| 74 | #include <iostream> |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 75 | #include <iterator> |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 76 | #include <fstream> |
| 77 | #include "ceres/ceres.h" |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 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"); |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 84 | DEFINE_string(minimizer, "trust_region", |
| 85 | "Minimizer type to use, choices are: line_search & trust_region"); |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 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"); |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 95 | DEFINE_string(line_search, "armijo", |
| 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."); |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 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 | |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 119 | namespace ceres { |
| 120 | namespace examples { |
| 121 | |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 122 | using Eigen::Dynamic; |
| 123 | using Eigen::RowMajor; |
| 124 | typedef Eigen::Matrix<double, Dynamic, 1> Vector; |
| 125 | typedef Eigen::Matrix<double, Dynamic, Dynamic, RowMajor> Matrix; |
| 126 | |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 127 | void 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 | |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 147 | bool GetAndSplitLine(std::ifstream& ifs, std::vector<std::string>* pieces) { |
| 148 | pieces->clear(); |
| 149 | char buf[256]; |
| 150 | ifs.getline(buf, 256); |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 151 | SplitStringUsingChar(std::string(buf), ' ', pieces); |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 152 | return true; |
| 153 | } |
| 154 | |
| 155 | void 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 | |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 162 | class NISTProblem { |
| 163 | public: |
| 164 | explicit NISTProblem(const std::string& filename) { |
| 165 | std::ifstream ifs(filename.c_str(), std::ifstream::in); |
| 166 | |
| 167 | std::vector<std::string> pieces; |
| 168 | SkipLines(ifs, 24); |
| 169 | GetAndSplitLine(ifs, &pieces); |
| 170 | const int kNumResponses = std::atoi(pieces[1].c_str()); |
| 171 | |
| 172 | GetAndSplitLine(ifs, &pieces); |
| 173 | const int kNumPredictors = std::atoi(pieces[0].c_str()); |
| 174 | |
| 175 | GetAndSplitLine(ifs, &pieces); |
| 176 | const int kNumObservations = std::atoi(pieces[0].c_str()); |
| 177 | |
| 178 | SkipLines(ifs, 4); |
| 179 | GetAndSplitLine(ifs, &pieces); |
| 180 | const int kNumParameters = std::atoi(pieces[0].c_str()); |
| 181 | SkipLines(ifs, 8); |
| 182 | |
| 183 | // Get the first line of initial and final parameter values to |
| 184 | // determine the number of tries. |
| 185 | GetAndSplitLine(ifs, &pieces); |
| 186 | const int kNumTries = pieces.size() - 4; |
| 187 | |
| 188 | predictor_.resize(kNumObservations, kNumPredictors); |
| 189 | response_.resize(kNumObservations, kNumResponses); |
| 190 | initial_parameters_.resize(kNumTries, kNumParameters); |
| 191 | final_parameters_.resize(1, kNumParameters); |
| 192 | |
| 193 | // Parse the line for parameter b1. |
| 194 | int parameter_id = 0; |
| 195 | for (int i = 0; i < kNumTries; ++i) { |
| 196 | initial_parameters_(i, parameter_id) = std::atof(pieces[i + 2].c_str()); |
| 197 | } |
| 198 | final_parameters_(0, parameter_id) = std::atof(pieces[2 + kNumTries].c_str()); |
| 199 | |
| 200 | // Parse the remaining parameter lines. |
| 201 | for (int parameter_id = 1; parameter_id < kNumParameters; ++parameter_id) { |
| 202 | GetAndSplitLine(ifs, &pieces); |
| 203 | // b2, b3, .... |
| 204 | for (int i = 0; i < kNumTries; ++i) { |
| 205 | initial_parameters_(i, parameter_id) = std::atof(pieces[i + 2].c_str()); |
| 206 | } |
| 207 | final_parameters_(0, parameter_id) = std::atof(pieces[2 + kNumTries].c_str()); |
| 208 | } |
| 209 | |
| 210 | // Certfied cost |
| 211 | SkipLines(ifs, 1); |
| 212 | GetAndSplitLine(ifs, &pieces); |
| 213 | certified_cost_ = std::atof(pieces[4].c_str()) / 2.0; |
| 214 | |
| 215 | // Read the observations. |
| 216 | SkipLines(ifs, 18 - kNumParameters); |
| 217 | for (int i = 0; i < kNumObservations; ++i) { |
| 218 | GetAndSplitLine(ifs, &pieces); |
| 219 | // Response. |
| 220 | for (int j = 0; j < kNumResponses; ++j) { |
| 221 | response_(i, j) = std::atof(pieces[j].c_str()); |
| 222 | } |
| 223 | |
| 224 | // Predictor variables. |
| 225 | for (int j = 0; j < kNumPredictors; ++j) { |
| 226 | predictor_(i, j) = std::atof(pieces[j + kNumResponses].c_str()); |
| 227 | } |
| 228 | } |
| 229 | } |
| 230 | |
| 231 | Matrix initial_parameters(int start) const { return initial_parameters_.row(start); } |
| 232 | Matrix final_parameters() const { return final_parameters_; } |
| 233 | Matrix predictor() const { return predictor_; } |
| 234 | Matrix response() const { return response_; } |
| 235 | int predictor_size() const { return predictor_.cols(); } |
| 236 | int num_observations() const { return predictor_.rows(); } |
| 237 | int response_size() const { return response_.cols(); } |
| 238 | int num_parameters() const { return initial_parameters_.cols(); } |
| 239 | int num_starts() const { return initial_parameters_.rows(); } |
| 240 | double certified_cost() const { return certified_cost_; } |
| 241 | |
| 242 | private: |
| 243 | Matrix predictor_; |
| 244 | Matrix response_; |
| 245 | Matrix initial_parameters_; |
| 246 | Matrix final_parameters_; |
| 247 | double certified_cost_; |
| 248 | }; |
| 249 | |
| 250 | #define NIST_BEGIN(CostFunctionName) \ |
| 251 | struct CostFunctionName { \ |
| 252 | CostFunctionName(const double* const x, \ |
| 253 | const double* const y) \ |
| 254 | : x_(*x), y_(*y) {} \ |
| 255 | double x_; \ |
| 256 | double y_; \ |
| 257 | template <typename T> \ |
| 258 | bool operator()(const T* const b, T* residual) const { \ |
| 259 | const T y(y_); \ |
| 260 | const T x(x_); \ |
| 261 | residual[0] = y - ( |
| 262 | |
| 263 | #define NIST_END ); return true; }}; |
| 264 | |
| 265 | // y = b1 * (b2+x)**(-1/b3) + e |
| 266 | NIST_BEGIN(Bennet5) |
| 267 | b[0] * pow(b[1] + x, T(-1.0) / b[2]) |
| 268 | NIST_END |
| 269 | |
| 270 | // y = b1*(1-exp[-b2*x]) + e |
| 271 | NIST_BEGIN(BoxBOD) |
| 272 | b[0] * (T(1.0) - exp(-b[1] * x)) |
| 273 | NIST_END |
| 274 | |
| 275 | // y = exp[-b1*x]/(b2+b3*x) + e |
| 276 | NIST_BEGIN(Chwirut) |
| 277 | exp(-b[0] * x) / (b[1] + b[2] * x) |
| 278 | NIST_END |
| 279 | |
| 280 | // y = b1*x**b2 + e |
| 281 | NIST_BEGIN(DanWood) |
| 282 | b[0] * pow(x, b[1]) |
| 283 | NIST_END |
| 284 | |
| 285 | // y = b1*exp( -b2*x ) + b3*exp( -(x-b4)**2 / b5**2 ) |
| 286 | // + b6*exp( -(x-b7)**2 / b8**2 ) + e |
| 287 | NIST_BEGIN(Gauss) |
| 288 | b[0] * exp(-b[1] * x) + |
| 289 | b[2] * exp(-pow((x - b[3])/b[4], 2)) + |
| 290 | b[5] * exp(-pow((x - b[6])/b[7],2)) |
| 291 | NIST_END |
| 292 | |
| 293 | // y = b1*exp(-b2*x) + b3*exp(-b4*x) + b5*exp(-b6*x) + e |
| 294 | NIST_BEGIN(Lanczos) |
| 295 | b[0] * exp(-b[1] * x) + b[2] * exp(-b[3] * x) + b[4] * exp(-b[5] * x) |
| 296 | NIST_END |
| 297 | |
| 298 | // y = (b1+b2*x+b3*x**2+b4*x**3) / |
| 299 | // (1+b5*x+b6*x**2+b7*x**3) + e |
| 300 | NIST_BEGIN(Hahn1) |
| 301 | (b[0] + b[1] * x + b[2] * x * x + b[3] * x * x * x) / |
| 302 | (T(1.0) + b[4] * x + b[5] * x * x + b[6] * x * x * x) |
| 303 | NIST_END |
| 304 | |
| 305 | // y = (b1 + b2*x + b3*x**2) / |
| 306 | // (1 + b4*x + b5*x**2) + e |
| 307 | NIST_BEGIN(Kirby2) |
| 308 | (b[0] + b[1] * x + b[2] * x * x) / |
| 309 | (T(1.0) + b[3] * x + b[4] * x * x) |
| 310 | NIST_END |
| 311 | |
| 312 | // y = b1*(x**2+x*b2) / (x**2+x*b3+b4) + e |
| 313 | NIST_BEGIN(MGH09) |
| 314 | b[0] * (x * x + x * b[1]) / (x * x + x * b[2] + b[3]) |
| 315 | NIST_END |
| 316 | |
| 317 | // y = b1 * exp[b2/(x+b3)] + e |
| 318 | NIST_BEGIN(MGH10) |
| 319 | b[0] * exp(b[1] / (x + b[2])) |
| 320 | NIST_END |
| 321 | |
| 322 | // y = b1 + b2*exp[-x*b4] + b3*exp[-x*b5] |
| 323 | NIST_BEGIN(MGH17) |
| 324 | b[0] + b[1] * exp(-x * b[3]) + b[2] * exp(-x * b[4]) |
| 325 | NIST_END |
| 326 | |
| 327 | // y = b1*(1-exp[-b2*x]) + e |
| 328 | NIST_BEGIN(Misra1a) |
| 329 | b[0] * (T(1.0) - exp(-b[1] * x)) |
| 330 | NIST_END |
| 331 | |
| 332 | // y = b1 * (1-(1+b2*x/2)**(-2)) + e |
| 333 | NIST_BEGIN(Misra1b) |
| 334 | b[0] * (T(1.0) - T(1.0)/ ((T(1.0) + b[1] * x / 2.0) * (T(1.0) + b[1] * x / 2.0))) |
| 335 | NIST_END |
| 336 | |
| 337 | // y = b1 * (1-(1+2*b2*x)**(-.5)) + e |
| 338 | NIST_BEGIN(Misra1c) |
| 339 | b[0] * (T(1.0) - pow(T(1.0) + T(2.0) * b[1] * x, -0.5)) |
| 340 | NIST_END |
| 341 | |
| 342 | // y = b1*b2*x*((1+b2*x)**(-1)) + e |
| 343 | NIST_BEGIN(Misra1d) |
| 344 | b[0] * b[1] * x / (T(1.0) + b[1] * x) |
| 345 | NIST_END |
| 346 | |
| 347 | const double kPi = 3.141592653589793238462643383279; |
| 348 | // pi = 3.141592653589793238462643383279E0 |
| 349 | // y = b1 - b2*x - arctan[b3/(x-b4)]/pi + e |
| 350 | NIST_BEGIN(Roszman1) |
| 351 | b[0] - b[1] * x - atan2(b[2], (x - b[3]))/T(kPi) |
| 352 | NIST_END |
| 353 | |
| 354 | // y = b1 / (1+exp[b2-b3*x]) + e |
| 355 | NIST_BEGIN(Rat42) |
| 356 | b[0] / (T(1.0) + exp(b[1] - b[2] * x)) |
| 357 | NIST_END |
| 358 | |
| 359 | // y = b1 / ((1+exp[b2-b3*x])**(1/b4)) + e |
| 360 | NIST_BEGIN(Rat43) |
| 361 | b[0] / pow(T(1.0) + exp(b[1] - b[2] * x), T(1.0) / b[3]) |
| 362 | NIST_END |
| 363 | |
| 364 | // y = (b1 + b2*x + b3*x**2 + b4*x**3) / |
| 365 | // (1 + b5*x + b6*x**2 + b7*x**3) + e |
| 366 | NIST_BEGIN(Thurber) |
| 367 | (b[0] + b[1] * x + b[2] * x * x + b[3] * x * x * x) / |
| 368 | (T(1.0) + b[4] * x + b[5] * x * x + b[6] * x * x * x) |
| 369 | NIST_END |
| 370 | |
| 371 | // y = b1 + b2*cos( 2*pi*x/12 ) + b3*sin( 2*pi*x/12 ) |
| 372 | // + b5*cos( 2*pi*x/b4 ) + b6*sin( 2*pi*x/b4 ) |
| 373 | // + b8*cos( 2*pi*x/b7 ) + b9*sin( 2*pi*x/b7 ) + e |
| 374 | NIST_BEGIN(ENSO) |
| 375 | b[0] + b[1] * cos(T(2.0 * kPi) * x / T(12.0)) + |
| 376 | b[2] * sin(T(2.0 * kPi) * x / T(12.0)) + |
| 377 | b[4] * cos(T(2.0 * kPi) * x / b[3]) + |
| 378 | b[5] * sin(T(2.0 * kPi) * x / b[3]) + |
| 379 | b[7] * cos(T(2.0 * kPi) * x / b[6]) + |
| 380 | b[8] * sin(T(2.0 * kPi) * x / b[6]) |
| 381 | NIST_END |
| 382 | |
| 383 | // y = (b1/b2) * exp[-0.5*((x-b3)/b2)**2] + e |
| 384 | NIST_BEGIN(Eckerle4) |
| 385 | b[0] / b[1] * exp(T(-0.5) * pow((x - b[2])/b[1], 2)) |
| 386 | NIST_END |
| 387 | |
| 388 | struct Nelson { |
| 389 | public: |
| 390 | Nelson(const double* const x, const double* const y) |
| 391 | : x1_(x[0]), x2_(x[1]), y_(y[0]) {} |
| 392 | |
| 393 | template <typename T> |
| 394 | bool operator()(const T* const b, T* residual) const { |
| 395 | // log[y] = b1 - b2*x1 * exp[-b3*x2] + e |
| 396 | residual[0] = T(log(y_)) - (b[0] - b[1] * T(x1_) * exp(-b[2] * T(x2_))); |
| 397 | return true; |
| 398 | } |
| 399 | |
| 400 | private: |
| 401 | double x1_; |
| 402 | double x2_; |
| 403 | double y_; |
| 404 | }; |
| 405 | |
| 406 | template <typename Model, int num_residuals, int num_parameters> |
| 407 | int RegressionDriver(const std::string& filename, |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 408 | const ceres::Solver::Options& options) { |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 409 | NISTProblem nist_problem(FLAGS_nist_data_dir + filename); |
| 410 | CHECK_EQ(num_residuals, nist_problem.response_size()); |
| 411 | CHECK_EQ(num_parameters, nist_problem.num_parameters()); |
| 412 | |
| 413 | Matrix predictor = nist_problem.predictor(); |
| 414 | Matrix response = nist_problem.response(); |
| 415 | Matrix final_parameters = nist_problem.final_parameters(); |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 416 | |
| 417 | printf("%s\n", filename.c_str()); |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 418 | |
| 419 | // Each NIST problem comes with multiple starting points, so we |
| 420 | // construct the problem from scratch for each case and solve it. |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 421 | int num_success = 0; |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 422 | for (int start = 0; start < nist_problem.num_starts(); ++start) { |
| 423 | Matrix initial_parameters = nist_problem.initial_parameters(start); |
| 424 | |
| 425 | ceres::Problem problem; |
| 426 | for (int i = 0; i < nist_problem.num_observations(); ++i) { |
| 427 | problem.AddResidualBlock( |
| 428 | new ceres::AutoDiffCostFunction<Model, num_residuals, num_parameters>( |
| 429 | new Model(predictor.data() + nist_problem.predictor_size() * i, |
| 430 | response.data() + nist_problem.response_size() * i)), |
| 431 | NULL, |
| 432 | initial_parameters.data()); |
| 433 | } |
| 434 | |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 435 | ceres::Solver::Summary summary; |
| 436 | Solve(options, &problem, &summary); |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 437 | |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 438 | // Compute the LRE by comparing each component of the solution |
| 439 | // with the ground truth, and taking the minimum. |
| 440 | Matrix final_parameters = nist_problem.final_parameters(); |
| 441 | const double kMaxNumSignificantDigits = 11; |
| 442 | double log_relative_error = kMaxNumSignificantDigits + 1; |
| 443 | for (int i = 0; i < num_parameters; ++i) { |
| 444 | const double tmp_lre = |
| 445 | -std::log10(std::fabs(final_parameters(i) - initial_parameters(i)) / |
| 446 | std::fabs(final_parameters(i))); |
| 447 | // The maximum LRE is capped at 11 - the precision at which the |
| 448 | // ground truth is known. |
| 449 | // |
| 450 | // The minimum LRE is capped at 0 - no digits match between the |
| 451 | // computed solution and the ground truth. |
| 452 | log_relative_error = |
| 453 | std::min(log_relative_error, |
| 454 | std::max(0.0, std::min(kMaxNumSignificantDigits, tmp_lre))); |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 455 | } |
| 456 | |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 457 | const int kMinNumMatchingDigits = 4; |
| 458 | if (log_relative_error >= kMinNumMatchingDigits) { |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 459 | ++num_success; |
| 460 | } |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 461 | |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 462 | printf("start: %d status: %s lre: %4.1f initial cost: %e final cost:%e " |
| 463 | "certified cost: %e total iterations: %d\n", |
| 464 | start + 1, |
| 465 | log_relative_error < kMinNumMatchingDigits ? "FAILURE" : "SUCCESS", |
| 466 | log_relative_error, |
| 467 | summary.initial_cost, |
| 468 | summary.final_cost, |
| 469 | nist_problem.certified_cost(), |
| 470 | (summary.num_successful_steps + summary.num_unsuccessful_steps)); |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 471 | } |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 472 | return num_success; |
| 473 | } |
| 474 | |
| 475 | void SetMinimizerOptions(ceres::Solver::Options* options) { |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 476 | CHECK(ceres::StringToMinimizerType(FLAGS_minimizer, |
| 477 | &options->minimizer_type)); |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 478 | CHECK(ceres::StringToLinearSolverType(FLAGS_linear_solver, |
| 479 | &options->linear_solver_type)); |
| 480 | CHECK(ceres::StringToPreconditionerType(FLAGS_preconditioner, |
| 481 | &options->preconditioner_type)); |
| 482 | CHECK(ceres::StringToTrustRegionStrategyType( |
| 483 | FLAGS_trust_region_strategy, |
| 484 | &options->trust_region_strategy_type)); |
| 485 | CHECK(ceres::StringToDoglegType(FLAGS_dogleg, &options->dogleg_type)); |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 486 | CHECK(ceres::StringToLineSearchDirectionType( |
| 487 | FLAGS_line_search_direction, |
| 488 | &options->line_search_direction_type)); |
| 489 | CHECK(ceres::StringToLineSearchType(FLAGS_line_search, |
| 490 | &options->line_search_type)); |
| 491 | CHECK(ceres::StringToLineSearchInterpolationType( |
| 492 | FLAGS_line_search_interpolation, |
| 493 | &options->line_search_interpolation_type)); |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 494 | |
| 495 | options->max_num_iterations = FLAGS_num_iterations; |
| 496 | options->use_nonmonotonic_steps = FLAGS_nonmonotonic_steps; |
| 497 | options->initial_trust_region_radius = FLAGS_initial_trust_region_radius; |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 498 | options->max_lbfgs_rank = FLAGS_lbfgs_rank; |
| 499 | options->line_search_sufficient_function_decrease = FLAGS_sufficient_decrease; |
| 500 | options->line_search_sufficient_curvature_decrease = |
| 501 | FLAGS_sufficient_curvature_decrease; |
| 502 | options->max_num_line_search_step_size_iterations = |
| 503 | FLAGS_max_line_search_iterations; |
| 504 | options->max_num_line_search_direction_restarts = |
| 505 | FLAGS_max_line_search_restarts; |
| 506 | options->use_approximate_eigenvalue_bfgs_scaling = |
| 507 | FLAGS_approximate_eigenvalue_bfgs_scaling; |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 508 | options->function_tolerance = 1e-18; |
| 509 | options->gradient_tolerance = 1e-18; |
| 510 | options->parameter_tolerance = 1e-18; |
| 511 | } |
| 512 | |
| 513 | void SolveNISTProblems() { |
| 514 | if (FLAGS_nist_data_dir.empty()) { |
| 515 | LOG(FATAL) << "Must specify the directory containing the NIST problems"; |
| 516 | } |
| 517 | |
| 518 | ceres::Solver::Options options; |
| 519 | SetMinimizerOptions(&options); |
| 520 | |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 521 | std::cout << "Lower Difficulty\n"; |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 522 | int easy_success = 0; |
| 523 | easy_success += RegressionDriver<Misra1a, 1, 2>("Misra1a.dat", options); |
| 524 | easy_success += RegressionDriver<Chwirut, 1, 3>("Chwirut1.dat", options); |
| 525 | easy_success += RegressionDriver<Chwirut, 1, 3>("Chwirut2.dat", options); |
| 526 | easy_success += RegressionDriver<Lanczos, 1, 6>("Lanczos3.dat", options); |
| 527 | easy_success += RegressionDriver<Gauss, 1, 8>("Gauss1.dat", options); |
| 528 | easy_success += RegressionDriver<Gauss, 1, 8>("Gauss2.dat", options); |
| 529 | easy_success += RegressionDriver<DanWood, 1, 2>("DanWood.dat", options); |
| 530 | easy_success += RegressionDriver<Misra1b, 1, 2>("Misra1b.dat", options); |
| 531 | |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 532 | std::cout << "\nMedium Difficulty\n"; |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 533 | int medium_success = 0; |
| 534 | medium_success += RegressionDriver<Kirby2, 1, 5>("Kirby2.dat", options); |
| 535 | medium_success += RegressionDriver<Hahn1, 1, 7>("Hahn1.dat", options); |
| 536 | medium_success += RegressionDriver<Nelson, 1, 3>("Nelson.dat", options); |
| 537 | medium_success += RegressionDriver<MGH17, 1, 5>("MGH17.dat", options); |
| 538 | medium_success += RegressionDriver<Lanczos, 1, 6>("Lanczos1.dat", options); |
| 539 | medium_success += RegressionDriver<Lanczos, 1, 6>("Lanczos2.dat", options); |
| 540 | medium_success += RegressionDriver<Gauss, 1, 8>("Gauss3.dat", options); |
| 541 | medium_success += RegressionDriver<Misra1c, 1, 2>("Misra1c.dat", options); |
| 542 | medium_success += RegressionDriver<Misra1d, 1, 2>("Misra1d.dat", options); |
| 543 | medium_success += RegressionDriver<Roszman1, 1, 4>("Roszman1.dat", options); |
| 544 | medium_success += RegressionDriver<ENSO, 1, 9>("ENSO.dat", options); |
| 545 | |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 546 | std::cout << "\nHigher Difficulty\n"; |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 547 | int hard_success = 0; |
| 548 | hard_success += RegressionDriver<MGH09, 1, 4>("MGH09.dat", options); |
| 549 | hard_success += RegressionDriver<Thurber, 1, 7>("Thurber.dat", options); |
| 550 | hard_success += RegressionDriver<BoxBOD, 1, 2>("BoxBOD.dat", options); |
| 551 | hard_success += RegressionDriver<Rat42, 1, 3>("Rat42.dat", options); |
| 552 | hard_success += RegressionDriver<MGH10, 1, 3>("MGH10.dat", options); |
| 553 | |
| 554 | hard_success += RegressionDriver<Eckerle4, 1, 3>("Eckerle4.dat", options); |
| 555 | hard_success += RegressionDriver<Rat43, 1, 4>("Rat43.dat", options); |
| 556 | hard_success += RegressionDriver<Bennet5, 1, 3>("Bennett5.dat", options); |
| 557 | |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 558 | std::cout << "\n"; |
| 559 | std::cout << "Easy : " << easy_success << "/16\n"; |
| 560 | std::cout << "Medium : " << medium_success << "/22\n"; |
| 561 | std::cout << "Hard : " << hard_success << "/16\n"; |
| 562 | std::cout << "Total : " << easy_success + medium_success + hard_success << "/54\n"; |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 563 | } |
| 564 | |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 565 | } // namespace examples |
| 566 | } // namespace ceres |
| 567 | |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 568 | int main(int argc, char** argv) { |
| 569 | google::ParseCommandLineFlags(&argc, &argv, true); |
| 570 | google::InitGoogleLogging(argv[0]); |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 571 | ceres::examples::SolveNISTProblems(); |
Angus Kong | 0ae28bd | 2013-02-13 14:56:04 -0800 | [diff] [blame] | 572 | return 0; |
| 573 | }; |