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