internal/ceres/solver.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: keir@google.com (Keir Mierle)
  30 //         sameeragarwal@google.com (Sameer Agarwal)
  31
  32 #include "ceres/solver.h"
  33
  34 #include <algorithm>
  35 #include <sstream>   // NOLINT
  36 #include <vector>
  37 #include "ceres/detect_structure.h"
  38 #include "ceres/gradient_checking_cost_function.h"
  39 #include "ceres/internal/port.h"
  40 #include "ceres/parameter_block_ordering.h"
  41 #include "ceres/preprocessor.h"
  42 #include "ceres/problem.h"
  43 #include "ceres/problem_impl.h"
  44 #include "ceres/program.h"
  45 #include "ceres/schur_templates.h"
  46 #include "ceres/solver_utils.h"
  47 #include "ceres/stringprintf.h"
  48 #include "ceres/types.h"
  49 #include "ceres/wall_time.h"
  50
  51 namespace ceres {
  52 namespace {
  53
  54 using std::map;
  55 using std::string;
  56 using std::vector;
  57
  58 #define OPTION_OP(x, y, OP)                                             \
  59   if (!(options.x OP y)) {                                              \
  60     std::stringstream ss;                                               \
  61     ss << "Invalid configuration. ";                                    \
  62     ss << string("Solver::Options::" #x " = ") << options.x << ". ";    \
  63     ss << "Violated constraint: ";                                      \
  64     ss << string("Solver::Options::" #x " " #OP " "#y);                 \
  65     *error = ss.str();                                                  \
  66     return false;                                                       \
  67   }
  68
  69 #define OPTION_OP_OPTION(x, y, OP)                                      \
  70   if (!(options.x OP options.y)) {                                      \
  71     std::stringstream ss;                                               \
  72     ss << "Invalid configuration. ";                                    \
  73     ss << string("Solver::Options::" #x " = ") << options.x << ". ";    \
  74     ss << string("Solver::Options::" #y " = ") << options.y << ". ";    \
  75     ss << "Violated constraint: ";                                      \
  76     ss << string("Solver::Options::" #x);                               \
  77     ss << string(#OP " Solver::Options::" #y ".");                      \
  78     *error = ss.str();                                                  \
  79     return false;                                                       \
  80   }
  81
  82 #define OPTION_GE(x, y) OPTION_OP(x, y, >=);
  83 #define OPTION_GT(x, y) OPTION_OP(x, y, >);
  84 #define OPTION_LE(x, y) OPTION_OP(x, y, <=);
  85 #define OPTION_LT(x, y) OPTION_OP(x, y, <);
  86 #define OPTION_LE_OPTION(x, y) OPTION_OP_OPTION(x, y, <=)
  87 #define OPTION_LT_OPTION(x, y) OPTION_OP_OPTION(x, y, <)
  88
  89 bool CommonOptionsAreValid(const Solver::Options& options, string* error) {
  90   OPTION_GE(max_num_iterations, 0);
  91   OPTION_GE(max_solver_time_in_seconds, 0.0);
  92   OPTION_GE(function_tolerance, 0.0);
  93   OPTION_GE(gradient_tolerance, 0.0);
  94   OPTION_GE(parameter_tolerance, 0.0);
  95   OPTION_GT(num_threads, 0);
  96   OPTION_GT(num_linear_solver_threads, 0);
  97   if (options.check_gradients) {
  98     OPTION_GT(gradient_check_relative_precision, 0.0);
  99     OPTION_GT(gradient_check_numeric_derivative_relative_step_size, 0.0);
 100   }
 101   return true;
 102 }
 103
 104 bool TrustRegionOptionsAreValid(const Solver::Options& options, string* error) {
 105   OPTION_GT(initial_trust_region_radius, 0.0);
 106   OPTION_GT(min_trust_region_radius, 0.0);
 107   OPTION_GT(max_trust_region_radius, 0.0);
 108   OPTION_LE_OPTION(min_trust_region_radius, max_trust_region_radius);
 109   OPTION_LE_OPTION(min_trust_region_radius, initial_trust_region_radius);
 110   OPTION_LE_OPTION(initial_trust_region_radius, max_trust_region_radius);
 111   OPTION_GE(min_relative_decrease, 0.0);
 112   OPTION_GE(min_lm_diagonal, 0.0);
 113   OPTION_GE(max_lm_diagonal, 0.0);
 114   OPTION_LE_OPTION(min_lm_diagonal, max_lm_diagonal);
 115   OPTION_GE(max_num_consecutive_invalid_steps, 0);
 116   OPTION_GT(eta, 0.0);
 117   OPTION_GE(min_linear_solver_iterations, 0);
 118   OPTION_GE(max_linear_solver_iterations, 1);
 119   OPTION_LE_OPTION(min_linear_solver_iterations, max_linear_solver_iterations);
 120
 121   if (options.use_inner_iterations) {
 122     OPTION_GE(inner_iteration_tolerance, 0.0);
 123   }
 124
 125   if (options.use_nonmonotonic_steps) {
 126     OPTION_GT(max_consecutive_nonmonotonic_steps, 0);
 127   }
 128
 129   if (options.linear_solver_type == ITERATIVE_SCHUR &&
 130       options.use_explicit_schur_complement &&
 131       options.preconditioner_type != SCHUR_JACOBI) {
 132     *error =  "use_explicit_schur_complement only supports "
 133         "SCHUR_JACOBI as the preconditioner.";
 134     return false;
 135   }
 136
 137   if (options.preconditioner_type == CLUSTER_JACOBI &&
 138       options.sparse_linear_algebra_library_type != SUITE_SPARSE) {
 139     *error =  "CLUSTER_JACOBI requires "
 140         "Solver::Options::sparse_linear_algebra_library_type to be "
 141         "SUITE_SPARSE";
 142     return false;
 143   }
 144
 145   if (options.preconditioner_type == CLUSTER_TRIDIAGONAL &&
 146       options.sparse_linear_algebra_library_type != SUITE_SPARSE) {
 147     *error =  "CLUSTER_TRIDIAGONAL requires "
 148         "Solver::Options::sparse_linear_algebra_library_type to be "
 149         "SUITE_SPARSE";
 150     return false;
 151   }
 152
 153 #ifdef CERES_NO_LAPACK
 154   if (options.dense_linear_algebra_library_type == LAPACK) {
 155     if (options.linear_solver_type == DENSE_NORMAL_CHOLESKY) {
 156       *error = "Can't use DENSE_NORMAL_CHOLESKY with LAPACK because "
 157           "LAPACK was not enabled when Ceres was built.";
 158       return false;
 159     } else if (options.linear_solver_type == DENSE_QR) {
 160       *error = "Can't use DENSE_QR with LAPACK because "
 161           "LAPACK was not enabled when Ceres was built.";
 162       return false;
 163     } else if (options.linear_solver_type == DENSE_SCHUR) {
 164       *error = "Can't use DENSE_SCHUR with LAPACK because "
 165           "LAPACK was not enabled when Ceres was built.";
 166       return false;
 167     }
 168   }
 169 #endif
 170
 171 #ifdef CERES_NO_SUITESPARSE
 172   if (options.sparse_linear_algebra_library_type == SUITE_SPARSE) {
 173     if (options.linear_solver_type == SPARSE_NORMAL_CHOLESKY) {
 174       *error = "Can't use SPARSE_NORMAL_CHOLESKY with SUITESPARSE because "
 175              "SuiteSparse was not enabled when Ceres was built.";
 176       return false;
 177     } else if (options.linear_solver_type == SPARSE_SCHUR) {
 178       *error = "Can't use SPARSE_SCHUR with SUITESPARSE because "
 179           "SuiteSparse was not enabled when Ceres was built.";
 180       return false;
 181     } else if (options.preconditioner_type == CLUSTER_JACOBI) {
 182       *error =  "CLUSTER_JACOBI preconditioner not supported. "
 183           "SuiteSparse was not enabled when Ceres was built.";
 184       return false;
 185     } else if (options.preconditioner_type == CLUSTER_TRIDIAGONAL) {
 186       *error =  "CLUSTER_TRIDIAGONAL preconditioner not supported. "
 187           "SuiteSparse was not enabled when Ceres was built.";
 188     return false;
 189     }
 190   }
 191 #endif
 192
 193 #ifdef CERES_NO_CXSPARSE
 194   if (options.sparse_linear_algebra_library_type == CX_SPARSE) {
 195     if (options.linear_solver_type == SPARSE_NORMAL_CHOLESKY) {
 196       *error = "Can't use SPARSE_NORMAL_CHOLESKY with CX_SPARSE because "
 197              "CXSparse was not enabled when Ceres was built.";
 198       return false;
 199     } else if (options.linear_solver_type == SPARSE_SCHUR) {
 200       *error = "Can't use SPARSE_SCHUR with CX_SPARSE because "
 201           "CXSparse was not enabled when Ceres was built.";
 202       return false;
 203     }
 204   }
 205 #endif
 206
 207 #ifndef CERES_USE_EIGEN_SPARSE
 208   if (options.sparse_linear_algebra_library_type == EIGEN_SPARSE) {
 209     if (options.linear_solver_type == SPARSE_NORMAL_CHOLESKY) {
 210       *error = "Can't use SPARSE_NORMAL_CHOLESKY with EIGEN_SPARSE because "
 211           "Eigen's sparse linear algebra was not enabled when Ceres was "
 212           "built.";
 213       return false;
 214     } else if (options.linear_solver_type == SPARSE_SCHUR) {
 215       *error = "Can't use SPARSE_SCHUR with EIGEN_SPARSE because "
 216           "Eigen's sparse linear algebra was not enabled when Ceres was "
 217           "built.";
 218       return false;
 219     }
 220   }
 221 #endif
 222
 223   if (options.sparse_linear_algebra_library_type == NO_SPARSE) {
 224     if (options.linear_solver_type == SPARSE_NORMAL_CHOLESKY) {
 225       *error = "Can't use SPARSE_NORMAL_CHOLESKY as "
 226           "sparse_linear_algebra_library_type is NO_SPARSE.";
 227       return false;
 228     } else if (options.linear_solver_type == SPARSE_SCHUR) {
 229       *error = "Can't use SPARSE_SCHUR as "
 230           "sparse_linear_algebra_library_type is NO_SPARSE.";
 231       return false;
 232     }
 233   }
 234
 235   if (options.trust_region_strategy_type == DOGLEG) {
 236     if (options.linear_solver_type == ITERATIVE_SCHUR ||
 237         options.linear_solver_type == CGNR) {
 238       *error = "DOGLEG only supports exact factorization based linear "
 239           "solvers. If you want to use an iterative solver please "
 240           "use LEVENBERG_MARQUARDT as the trust_region_strategy_type";
 241       return false;
 242     }
 243   }
 244
 245   if (options.trust_region_minimizer_iterations_to_dump.size() > 0 &&
 246       options.trust_region_problem_dump_format_type != CONSOLE &&
 247       options.trust_region_problem_dump_directory.empty()) {
 248     *error = "Solver::Options::trust_region_problem_dump_directory is empty.";
 249     return false;
 250   }
 251
 252   if (options.dynamic_sparsity &&
 253       options.linear_solver_type != SPARSE_NORMAL_CHOLESKY) {
 254     *error = "Dynamic sparsity is only supported with SPARSE_NORMAL_CHOLESKY.";
 255     return false;
 256   }
 257
 258   return true;
 259 }
 260
 261 bool LineSearchOptionsAreValid(const Solver::Options& options, string* error) {
 262   OPTION_GT(max_lbfgs_rank, 0);
 263   OPTION_GT(min_line_search_step_size, 0.0);
 264   OPTION_GT(max_line_search_step_contraction, 0.0);
 265   OPTION_LT(max_line_search_step_contraction, 1.0);
 266   OPTION_LT_OPTION(max_line_search_step_contraction,
 267                    min_line_search_step_contraction);
 268   OPTION_LE(min_line_search_step_contraction, 1.0);
 269   OPTION_GT(max_num_line_search_step_size_iterations, 0);
 270   OPTION_GT(line_search_sufficient_function_decrease, 0.0);
 271   OPTION_LT_OPTION(line_search_sufficient_function_decrease,
 272                    line_search_sufficient_curvature_decrease);
 273   OPTION_LT(line_search_sufficient_curvature_decrease, 1.0);
 274   OPTION_GT(max_line_search_step_expansion, 1.0);
 275
 276   if ((options.line_search_direction_type == ceres::BFGS ||
 277        options.line_search_direction_type == ceres::LBFGS) &&
 278       options.line_search_type != ceres::WOLFE) {
 279     *error =
 280         string("Invalid configuration: Solver::Options::line_search_type = ")
 281         + string(LineSearchTypeToString(options.line_search_type))
 282         + string(". When using (L)BFGS, "
 283                  "Solver::Options::line_search_type must be set to WOLFE.");
 284     return false;
 285   }
 286
 287   // Warn user if they have requested BISECTION interpolation, but constraints
 288   // on max/min step size change during line search prevent bisection scaling
 289   // from occurring. Warn only, as this is likely a user mistake, but one which
 290   // does not prevent us from continuing.
 291   LOG_IF(WARNING,
 292          (options.line_search_interpolation_type == ceres::BISECTION &&
 293           (options.max_line_search_step_contraction > 0.5 ||
 294            options.min_line_search_step_contraction < 0.5)))
 295       << "Line search interpolation type is BISECTION, but specified "
 296       << "max_line_search_step_contraction: "
 297       << options.max_line_search_step_contraction << ", and "
 298       << "min_line_search_step_contraction: "
 299       << options.min_line_search_step_contraction
 300       << ", prevent bisection (0.5) scaling, continuing with solve regardless.";
 301
 302   return true;
 303 }
 304
 305 #undef OPTION_OP
 306 #undef OPTION_OP_OPTION
 307 #undef OPTION_GT
 308 #undef OPTION_GE
 309 #undef OPTION_LE
 310 #undef OPTION_LT
 311 #undef OPTION_LE_OPTION
 312 #undef OPTION_LT_OPTION
 313
 314 void StringifyOrdering(const vector<int>& ordering, string* report) {
 315   if (ordering.size() == 0) {
 316     internal::StringAppendF(report, "AUTOMATIC");
 317     return;
 318   }
 319
 320   for (int i = 0; i < ordering.size() - 1; ++i) {
 321     internal::StringAppendF(report, "%d,", ordering[i]);
 322   }
 323   internal::StringAppendF(report, "%d", ordering.back());
 324 }
 325
 326 void SummarizeGivenProgram(const internal::Program& program,
 327                            Solver::Summary* summary) {
 328   summary->num_parameter_blocks     = program.NumParameterBlocks();
 329   summary->num_parameters           = program.NumParameters();
 330   summary->num_effective_parameters = program.NumEffectiveParameters();
 331   summary->num_residual_blocks      = program.NumResidualBlocks();
 332   summary->num_residuals            = program.NumResiduals();
 333 }
 334
 335 void SummarizeReducedProgram(const internal::Program& program,
 336                              Solver::Summary* summary) {
 337   summary->num_parameter_blocks_reduced     = program.NumParameterBlocks();
 338   summary->num_parameters_reduced           = program.NumParameters();
 339   summary->num_effective_parameters_reduced = program.NumEffectiveParameters();
 340   summary->num_residual_blocks_reduced      = program.NumResidualBlocks();
 341   summary->num_residuals_reduced            = program.NumResiduals();
 342 }
 343
 344 void PreSolveSummarize(const Solver::Options& options,
 345                        const internal::ProblemImpl* problem,
 346                        Solver::Summary* summary) {
 347   SummarizeGivenProgram(problem->program(), summary);
 348   internal::OrderingToGroupSizes(options.linear_solver_ordering.get(),
 349                                  &(summary->linear_solver_ordering_given));
 350   internal::OrderingToGroupSizes(options.inner_iteration_ordering.get(),
 351                                  &(summary->inner_iteration_ordering_given));
 352
 353   summary->dense_linear_algebra_library_type  = options.dense_linear_algebra_library_type;  //  NOLINT
 354   summary->dogleg_type                        = options.dogleg_type;
 355   summary->inner_iteration_time_in_seconds    = 0.0;
 356   summary->num_line_search_steps              = 0;
 357   summary->line_search_cost_evaluation_time_in_seconds = 0.0;
 358   summary->line_search_gradient_evaluation_time_in_seconds = 0.0;
 359   summary->line_search_polynomial_minimization_time_in_seconds = 0.0;
 360   summary->line_search_total_time_in_seconds  = 0.0;
 361   summary->inner_iterations_given             = options.use_inner_iterations;
 362   summary->line_search_direction_type         = options.line_search_direction_type;         //  NOLINT
 363   summary->line_search_interpolation_type     = options.line_search_interpolation_type;     //  NOLINT
 364   summary->line_search_type                   = options.line_search_type;
 365   summary->linear_solver_type_given           = options.linear_solver_type;
 366   summary->max_lbfgs_rank                     = options.max_lbfgs_rank;
 367   summary->minimizer_type                     = options.minimizer_type;
 368   summary->nonlinear_conjugate_gradient_type  = options.nonlinear_conjugate_gradient_type;  //  NOLINT
 369   summary->num_linear_solver_threads_given    = options.num_linear_solver_threads;          //  NOLINT
 370   summary->num_threads_given                  = options.num_threads;
 371   summary->preconditioner_type_given          = options.preconditioner_type;
 372   summary->sparse_linear_algebra_library_type = options.sparse_linear_algebra_library_type; //  NOLINT
 373   summary->trust_region_strategy_type         = options.trust_region_strategy_type;         //  NOLINT
 374   summary->visibility_clustering_type         = options.visibility_clustering_type;         //  NOLINT
 375 }
 376
 377 void PostSolveSummarize(const internal::PreprocessedProblem& pp,
 378                         Solver::Summary* summary) {
 379   internal::OrderingToGroupSizes(pp.options.linear_solver_ordering.get(),
 380                                  &(summary->linear_solver_ordering_used));
 381   internal::OrderingToGroupSizes(pp.options.inner_iteration_ordering.get(),
 382                                  &(summary->inner_iteration_ordering_used));
 383
 384   summary->inner_iterations_used          = pp.inner_iteration_minimizer.get() != NULL;     // NOLINT
 385   summary->linear_solver_type_used        = pp.linear_solver_options.type;
 386   summary->num_linear_solver_threads_used = pp.options.num_linear_solver_threads;           // NOLINT
 387   summary->num_threads_used               = pp.options.num_threads;
 388   summary->preconditioner_type_used       = pp.options.preconditioner_type;                 // NOLINT
 389
 390   internal::SetSummaryFinalCost(summary);
 391
 392   if (pp.reduced_program.get() != NULL) {
 393     SummarizeReducedProgram(*pp.reduced_program, summary);
 394   }
 395
 396   // It is possible that no evaluator was created. This would be the
 397   // case if the preprocessor failed, or if the reduced problem did
 398   // not contain any parameter blocks. Thus, only extract the
 399   // evaluator statistics if one exists.
 400   if (pp.evaluator.get() != NULL) {
 401     const map<string, double>& evaluator_time_statistics =
 402         pp.evaluator->TimeStatistics();
 403     summary->residual_evaluation_time_in_seconds =
 404         FindWithDefault(evaluator_time_statistics, "Evaluator::Residual", 0.0);
 405     summary->jacobian_evaluation_time_in_seconds =
 406         FindWithDefault(evaluator_time_statistics, "Evaluator::Jacobian", 0.0);
 407   }
 408
 409   // Again, like the evaluator, there may or may not be a linear
 410   // solver from which we can extract run time statistics. In
 411   // particular the line search solver does not use a linear solver.
 412   if (pp.linear_solver.get() != NULL) {
 413     const map<string, double>& linear_solver_time_statistics =
 414         pp.linear_solver->TimeStatistics();
 415     summary->linear_solver_time_in_seconds =
 416         FindWithDefault(linear_solver_time_statistics,
 417                         "LinearSolver::Solve",
 418                         0.0);
 419   }
 420 }
 421
 422 void Minimize(internal::PreprocessedProblem* pp,
 423               Solver::Summary* summary) {
 424   using internal::Program;
 425   using internal::scoped_ptr;
 426   using internal::Minimizer;
 427
 428   Program* program = pp->reduced_program.get();
 429   if (pp->reduced_program->NumParameterBlocks() == 0) {
 430     summary->message = "Function tolerance reached. "
 431         "No non-constant parameter blocks found.";
 432     summary->termination_type = CONVERGENCE;
 433     VLOG_IF(1, pp->options.logging_type != SILENT) << summary->message;
 434     summary->initial_cost = summary->fixed_cost;
 435     summary->final_cost = summary->fixed_cost;
 436     return;
 437   }
 438
 439   scoped_ptr<Minimizer> minimizer(
 440       Minimizer::Create(pp->options.minimizer_type));
 441   minimizer->Minimize(pp->minimizer_options,
 442                       pp->reduced_parameters.data(),
 443                       summary);
 444
 445   if (summary->IsSolutionUsable()) {
 446     program->StateVectorToParameterBlocks(pp->reduced_parameters.data());
 447     program->CopyParameterBlockStateToUserState();
 448   }
 449 }
 450
 451 std::string SchurStructureToString(const int row_block_size,
 452                                    const int e_block_size,
 453                                    const int f_block_size) {
 454   const std::string row =
 455       (row_block_size == Eigen::Dynamic)
 456       ? "d" : internal::StringPrintf("%d", row_block_size);
 457
 458   const std::string e =
 459       (e_block_size == Eigen::Dynamic)
 460       ? "d" : internal::StringPrintf("%d", e_block_size);
 461
 462   const std::string f =
 463       (f_block_size == Eigen::Dynamic)
 464       ? "d" : internal::StringPrintf("%d", f_block_size);
 465
 466   return internal::StringPrintf("%s,%s,%s", row.c_str(), e.c_str(), f.c_str());
 467 }
 468
 469 }  // namespace
 470
 471 bool Solver::Options::IsValid(string* error) const {
 472   if (!CommonOptionsAreValid(*this, error)) {
 473     return false;
 474   }
 475
 476   if (minimizer_type == TRUST_REGION &&
 477       !TrustRegionOptionsAreValid(*this, error)) {
 478     return false;
 479   }
 480
 481   // We do not know if the problem is bounds constrained or not, if it
 482   // is then the trust region solver will also use the line search
 483   // solver to do a projection onto the box constraints, so make sure
 484   // that the line search options are checked independent of what
 485   // minimizer algorithm is being used.
 486   return LineSearchOptionsAreValid(*this, error);
 487 }
 488
 489 Solver::~Solver() {}
 490
 491 void Solver::Solve(const Solver::Options& options,
 492                    Problem* problem,
 493                    Solver::Summary* summary) {
 494   using internal::PreprocessedProblem;
 495   using internal::Preprocessor;
 496   using internal::ProblemImpl;
 497   using internal::Program;
 498   using internal::scoped_ptr;
 499   using internal::WallTimeInSeconds;
 500
 501   CHECK_NOTNULL(problem);
 502   CHECK_NOTNULL(summary);
 503
 504   double start_time = WallTimeInSeconds();
 505   *summary = Summary();
 506   if (!options.IsValid(&summary->message)) {
 507     LOG(ERROR) << "Terminating: " << summary->message;
 508     return;
 509   }
 510
 511   ProblemImpl* problem_impl = problem->problem_impl_.get();
 512   Program* program = problem_impl->mutable_program();
 513   PreSolveSummarize(options, problem_impl, summary);
 514
 515   // Make sure that all the parameter blocks states are set to the
 516   // values provided by the user.
 517   program->SetParameterBlockStatePtrsToUserStatePtrs();
 518
 519   // If gradient_checking is enabled, wrap all cost functions in a
 520   // gradient checker and install a callback that terminates if any gradient
 521   // error is detected.
 522   scoped_ptr<internal::ProblemImpl> gradient_checking_problem;
 523   internal::GradientCheckingIterationCallback gradient_checking_callback;
 524   Solver::Options modified_options = options;
 525   if (options.check_gradients) {
 526     modified_options.callbacks.push_back(&gradient_checking_callback);
 527     gradient_checking_problem.reset(
 528         CreateGradientCheckingProblemImpl(
 529             problem_impl,
 530             options.gradient_check_numeric_derivative_relative_step_size,
 531             options.gradient_check_relative_precision,
 532             &gradient_checking_callback));
 533     problem_impl = gradient_checking_problem.get();
 534     program = problem_impl->mutable_program();
 535   }
 536
 537   scoped_ptr<Preprocessor> preprocessor(
 538       Preprocessor::Create(modified_options.minimizer_type));
 539   PreprocessedProblem pp;
 540
 541   const bool status = preprocessor->Preprocess(modified_options, problem_impl, &pp);
 542
 543   // We check the linear_solver_options.type rather than
 544   // modified_options.linear_solver_type because, depending on the
 545   // lack of a Schur structure, the preprocessor may change the linear
 546   // solver type.
 547   if (IsSchurType(pp.linear_solver_options.type)) {
 548     // TODO(sameeragarwal): We can likely eliminate the duplicate call
 549     // to DetectStructure here and inside the linear solver, by
 550     // calling this in the preprocessor.
 551     int row_block_size;
 552     int e_block_size;
 553     int f_block_size;
 554     DetectStructure(*static_cast<internal::BlockSparseMatrix*>(
 555                         pp.minimizer_options.jacobian.get())
 556                     ->block_structure(),
 557                     pp.linear_solver_options.elimination_groups[0],
 558                     &row_block_size,
 559                     &e_block_size,
 560                     &f_block_size);
 561     summary->schur_structure_given =
 562         SchurStructureToString(row_block_size, e_block_size, f_block_size);
 563     internal::GetBestSchurTemplateSpecialization(&row_block_size,
 564                                                  &e_block_size,
 565                                                  &f_block_size);
 566     summary->schur_structure_used =
 567         SchurStructureToString(row_block_size, e_block_size, f_block_size);
 568   }
 569
 570   summary->fixed_cost = pp.fixed_cost;
 571   summary->preprocessor_time_in_seconds = WallTimeInSeconds() - start_time;
 572
 573   if (status) {
 574     const double minimizer_start_time = WallTimeInSeconds();
 575     Minimize(&pp, summary);
 576     summary->minimizer_time_in_seconds =
 577         WallTimeInSeconds() - minimizer_start_time;
 578   } else {
 579     summary->message = pp.error;
 580   }
 581
 582   const double postprocessor_start_time = WallTimeInSeconds();
 583   problem_impl = problem->problem_impl_.get();
 584   program = problem_impl->mutable_program();
 585   // On exit, ensure that the parameter blocks again point at the user
 586   // provided values and the parameter blocks are numbered according
 587   // to their position in the original user provided program.
 588   program->SetParameterBlockStatePtrsToUserStatePtrs();
 589   program->SetParameterOffsetsAndIndex();
 590   PostSolveSummarize(pp, summary);
 591   summary->postprocessor_time_in_seconds =
 592       WallTimeInSeconds() - postprocessor_start_time;
 593
 594   // If the gradient checker reported an error, we want to report FAILURE
 595   // instead of USER_FAILURE and provide the error log.
 596   if (gradient_checking_callback.gradient_error_detected()) {
 597     summary->termination_type = FAILURE;
 598     summary->message = gradient_checking_callback.error_log();
 599   }
 600
 601   summary->total_time_in_seconds = WallTimeInSeconds() - start_time;
 602 }
 603
 604 void Solve(const Solver::Options& options,
 605            Problem* problem,
 606            Solver::Summary* summary) {
 607   Solver solver;
 608   solver.Solve(options, problem, summary);
 609 }
 610
 611 Solver::Summary::Summary()
 612     // Invalid values for most fields, to ensure that we are not
 613     // accidentally reporting default values.
 614     : minimizer_type(TRUST_REGION),
 615       termination_type(FAILURE),
 616       message("ceres::Solve was not called."),
 617       initial_cost(-1.0),
 618       final_cost(-1.0),
 619       fixed_cost(-1.0),
 620       num_successful_steps(-1),
 621       num_unsuccessful_steps(-1),
 622       num_inner_iteration_steps(-1),
 623       num_line_search_steps(-1),
 624       preprocessor_time_in_seconds(-1.0),
 625       minimizer_time_in_seconds(-1.0),
 626       postprocessor_time_in_seconds(-1.0),
 627       total_time_in_seconds(-1.0),
 628       linear_solver_time_in_seconds(-1.0),
 629       residual_evaluation_time_in_seconds(-1.0),
 630       jacobian_evaluation_time_in_seconds(-1.0),
 631       inner_iteration_time_in_seconds(-1.0),
 632       line_search_cost_evaluation_time_in_seconds(-1.0),
 633       line_search_gradient_evaluation_time_in_seconds(-1.0),
 634       line_search_polynomial_minimization_time_in_seconds(-1.0),
 635       line_search_total_time_in_seconds(-1.0),
 636       num_parameter_blocks(-1),
 637       num_parameters(-1),
 638       num_effective_parameters(-1),
 639       num_residual_blocks(-1),
 640       num_residuals(-1),
 641       num_parameter_blocks_reduced(-1),
 642       num_parameters_reduced(-1),
 643       num_effective_parameters_reduced(-1),
 644       num_residual_blocks_reduced(-1),
 645       num_residuals_reduced(-1),
 646       is_constrained(false),
 647       num_threads_given(-1),
 648       num_threads_used(-1),
 649       num_linear_solver_threads_given(-1),
 650       num_linear_solver_threads_used(-1),
 651       linear_solver_type_given(SPARSE_NORMAL_CHOLESKY),
 652       linear_solver_type_used(SPARSE_NORMAL_CHOLESKY),
 653       inner_iterations_given(false),
 654       inner_iterations_used(false),
 655       preconditioner_type_given(IDENTITY),
 656       preconditioner_type_used(IDENTITY),
 657       visibility_clustering_type(CANONICAL_VIEWS),
 658       trust_region_strategy_type(LEVENBERG_MARQUARDT),
 659       dense_linear_algebra_library_type(EIGEN),
 660       sparse_linear_algebra_library_type(SUITE_SPARSE),
 661       line_search_direction_type(LBFGS),
 662       line_search_type(ARMIJO),
 663       line_search_interpolation_type(BISECTION),
 664       nonlinear_conjugate_gradient_type(FLETCHER_REEVES),
 665       max_lbfgs_rank(-1) {
 666 }
 667
 668 using internal::StringAppendF;
 669 using internal::StringPrintf;
 670
 671 string Solver::Summary::BriefReport() const {
 672   return StringPrintf("Ceres Solver Report: "
 673                       "Iterations: %d, "
 674                       "Initial cost: %e, "
 675                       "Final cost: %e, "
 676                       "Termination: %s",
 677                       num_successful_steps + num_unsuccessful_steps,
 678                       initial_cost,
 679                       final_cost,
 680                       TerminationTypeToString(termination_type));
 681 }
 682
 683 string Solver::Summary::FullReport() const {
 684   using internal::VersionString;
 685
 686   string report = string("\nSolver Summary (v " + VersionString() + ")\n\n");
 687
 688   StringAppendF(&report, "%45s    %21s\n", "Original", "Reduced");
 689   StringAppendF(&report, "Parameter blocks    % 25d% 25d\n",
 690                 num_parameter_blocks, num_parameter_blocks_reduced);
 691   StringAppendF(&report, "Parameters          % 25d% 25d\n",
 692                 num_parameters, num_parameters_reduced);
 693   if (num_effective_parameters_reduced != num_parameters_reduced) {
 694     StringAppendF(&report, "Effective parameters% 25d% 25d\n",
 695                   num_effective_parameters, num_effective_parameters_reduced);
 696   }
 697   StringAppendF(&report, "Residual blocks     % 25d% 25d\n",
 698                 num_residual_blocks, num_residual_blocks_reduced);
 699   StringAppendF(&report, "Residual            % 25d% 25d\n",
 700                 num_residuals, num_residuals_reduced);
 701
 702   if (minimizer_type == TRUST_REGION) {
 703     // TRUST_SEARCH HEADER
 704     StringAppendF(&report, "\nMinimizer                 %19s\n",
 705                   "TRUST_REGION");
 706
 707     if (linear_solver_type_used == DENSE_NORMAL_CHOLESKY ||
 708         linear_solver_type_used == DENSE_SCHUR ||
 709         linear_solver_type_used == DENSE_QR) {
 710       StringAppendF(&report, "\nDense linear algebra library  %15s\n",
 711                     DenseLinearAlgebraLibraryTypeToString(
 712                         dense_linear_algebra_library_type));
 713     }
 714
 715     if (linear_solver_type_used == SPARSE_NORMAL_CHOLESKY ||
 716         linear_solver_type_used == SPARSE_SCHUR ||
 717         (linear_solver_type_used == ITERATIVE_SCHUR &&
 718          (preconditioner_type_used == CLUSTER_JACOBI ||
 719           preconditioner_type_used == CLUSTER_TRIDIAGONAL))) {
 720       StringAppendF(&report, "\nSparse linear algebra library %15s\n",
 721                     SparseLinearAlgebraLibraryTypeToString(
 722                         sparse_linear_algebra_library_type));
 723     }
 724
 725     StringAppendF(&report, "Trust region strategy     %19s",
 726                   TrustRegionStrategyTypeToString(
 727                       trust_region_strategy_type));
 728     if (trust_region_strategy_type == DOGLEG) {
 729       if (dogleg_type == TRADITIONAL_DOGLEG) {
 730         StringAppendF(&report, " (TRADITIONAL)");
 731       } else {
 732         StringAppendF(&report, " (SUBSPACE)");
 733       }
 734     }
 735     StringAppendF(&report, "\n");
 736     StringAppendF(&report, "\n");
 737
 738     StringAppendF(&report, "%45s    %21s\n", "Given",  "Used");
 739     StringAppendF(&report, "Linear solver       %25s%25s\n",
 740                   LinearSolverTypeToString(linear_solver_type_given),
 741                   LinearSolverTypeToString(linear_solver_type_used));
 742
 743     if (linear_solver_type_given == CGNR ||
 744         linear_solver_type_given == ITERATIVE_SCHUR) {
 745       StringAppendF(&report, "Preconditioner      %25s%25s\n",
 746                     PreconditionerTypeToString(preconditioner_type_given),
 747                     PreconditionerTypeToString(preconditioner_type_used));
 748     }
 749
 750     if (preconditioner_type_used == CLUSTER_JACOBI ||
 751         preconditioner_type_used == CLUSTER_TRIDIAGONAL) {
 752       StringAppendF(&report, "Visibility clustering%24s%25s\n",
 753                     VisibilityClusteringTypeToString(
 754                         visibility_clustering_type),
 755                     VisibilityClusteringTypeToString(
 756                         visibility_clustering_type));
 757     }
 758     StringAppendF(&report, "Threads             % 25d% 25d\n",
 759                   num_threads_given, num_threads_used);
 760     StringAppendF(&report, "Linear solver threads % 23d% 25d\n",
 761                   num_linear_solver_threads_given,
 762                   num_linear_solver_threads_used);
 763
 764     string given;
 765     StringifyOrdering(linear_solver_ordering_given, &given);
 766     string used;
 767     StringifyOrdering(linear_solver_ordering_used, &used);
 768     StringAppendF(&report,
 769                   "Linear solver ordering %22s %24s\n",
 770                   given.c_str(),
 771                   used.c_str());
 772     if (IsSchurType(linear_solver_type_used)) {
 773       StringAppendF(&report,
 774                     "Schur structure        %22s %24s\n",
 775                     schur_structure_given.c_str(),
 776                     schur_structure_used.c_str());
 777     }
 778
 779     if (inner_iterations_given) {
 780       StringAppendF(&report,
 781                     "Use inner iterations     %20s     %20s\n",
 782                     inner_iterations_given ? "True" : "False",
 783                     inner_iterations_used ? "True" : "False");
 784     }
 785
 786     if (inner_iterations_used) {
 787       string given;
 788       StringifyOrdering(inner_iteration_ordering_given, &given);
 789       string used;
 790       StringifyOrdering(inner_iteration_ordering_used, &used);
 791     StringAppendF(&report,
 792                   "Inner iteration ordering %20s %24s\n",
 793                   given.c_str(),
 794                   used.c_str());
 795     }
 796   } else {
 797     // LINE_SEARCH HEADER
 798     StringAppendF(&report, "\nMinimizer                 %19s\n", "LINE_SEARCH");
 799
 800
 801     string line_search_direction_string;
 802     if (line_search_direction_type == LBFGS) {
 803       line_search_direction_string = StringPrintf("LBFGS (%d)", max_lbfgs_rank);
 804     } else if (line_search_direction_type == NONLINEAR_CONJUGATE_GRADIENT) {
 805       line_search_direction_string =
 806           NonlinearConjugateGradientTypeToString(
 807               nonlinear_conjugate_gradient_type);
 808     } else {
 809       line_search_direction_string =
 810           LineSearchDirectionTypeToString(line_search_direction_type);
 811     }
 812
 813     StringAppendF(&report, "Line search direction     %19s\n",
 814                   line_search_direction_string.c_str());
 815
 816     const string line_search_type_string =
 817         StringPrintf("%s %s",
 818                      LineSearchInterpolationTypeToString(
 819                          line_search_interpolation_type),
 820                      LineSearchTypeToString(line_search_type));
 821     StringAppendF(&report, "Line search type          %19s\n",
 822                   line_search_type_string.c_str());
 823     StringAppendF(&report, "\n");
 824
 825     StringAppendF(&report, "%45s    %21s\n", "Given",  "Used");
 826     StringAppendF(&report, "Threads             % 25d% 25d\n",
 827                   num_threads_given, num_threads_used);
 828   }
 829
 830   StringAppendF(&report, "\nCost:\n");
 831   StringAppendF(&report, "Initial        % 30e\n", initial_cost);
 832   if (termination_type != FAILURE &&
 833       termination_type != USER_FAILURE) {
 834     StringAppendF(&report, "Final          % 30e\n", final_cost);
 835     StringAppendF(&report, "Change         % 30e\n",
 836                   initial_cost - final_cost);
 837   }
 838
 839   StringAppendF(&report, "\nMinimizer iterations         % 16d\n",
 840                 num_successful_steps + num_unsuccessful_steps);
 841
 842   // Successful/Unsuccessful steps only matter in the case of the
 843   // trust region solver. Line search terminates when it encounters
 844   // the first unsuccessful step.
 845   if (minimizer_type == TRUST_REGION) {
 846     StringAppendF(&report, "Successful steps               % 14d\n",
 847                   num_successful_steps);
 848     StringAppendF(&report, "Unsuccessful steps             % 14d\n",
 849                   num_unsuccessful_steps);
 850   }
 851   if (inner_iterations_used) {
 852     StringAppendF(&report, "Steps with inner iterations    % 14d\n",
 853                   num_inner_iteration_steps);
 854   }
 855
 856   const bool line_search_used =
 857       (minimizer_type == LINE_SEARCH ||
 858        (minimizer_type == TRUST_REGION && is_constrained));
 859
 860   if (line_search_used) {
 861     StringAppendF(&report, "Line search steps              % 14d\n",
 862                   num_line_search_steps);
 863   }
 864
 865   StringAppendF(&report, "\nTime (in seconds):\n");
 866   StringAppendF(&report, "Preprocessor        %25.4f\n",
 867                 preprocessor_time_in_seconds);
 868
 869   StringAppendF(&report, "\n  Residual evaluation %23.4f\n",
 870                 residual_evaluation_time_in_seconds);
 871   if (line_search_used) {
 872     StringAppendF(&report, "    Line search cost evaluation    %10.4f\n",
 873                   line_search_cost_evaluation_time_in_seconds);
 874   }
 875   StringAppendF(&report, "  Jacobian evaluation %23.4f\n",
 876                 jacobian_evaluation_time_in_seconds);
 877   if (line_search_used) {
 878     StringAppendF(&report, "    Line search gradient evaluation   %6.4f\n",
 879                   line_search_gradient_evaluation_time_in_seconds);
 880   }
 881
 882   if (minimizer_type == TRUST_REGION) {
 883     StringAppendF(&report, "  Linear solver       %23.4f\n",
 884                   linear_solver_time_in_seconds);
 885   }
 886
 887   if (inner_iterations_used) {
 888     StringAppendF(&report, "  Inner iterations    %23.4f\n",
 889                   inner_iteration_time_in_seconds);
 890   }
 891
 892   if (line_search_used) {
 893     StringAppendF(&report, "  Line search polynomial minimization  %.4f\n",
 894                   line_search_polynomial_minimization_time_in_seconds);
 895   }
 896
 897   StringAppendF(&report, "Minimizer           %25.4f\n\n",
 898                 minimizer_time_in_seconds);
 899
 900   StringAppendF(&report, "Postprocessor        %24.4f\n",
 901                 postprocessor_time_in_seconds);
 902
 903   StringAppendF(&report, "Total               %25.4f\n\n",
 904                 total_time_in_seconds);
 905
 906   StringAppendF(&report, "Termination:        %25s (%s)\n",
 907                 TerminationTypeToString(termination_type), message.c_str());
 908   return report;
 909 }
 910
 911 bool Solver::Summary::IsSolutionUsable() const {
 912   return internal::IsSolutionUsable(*this);
 913 }
 914
 915 }  // namespace ceres