9f3228bb0be19818ae5103befb7564fb9085ca55
[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(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->line_search_cost_evaluation_time_in_seconds = 0.0;
355   summary->line_search_gradient_evaluation_time_in_seconds = 0.0;
356   summary->line_search_polynomial_minimization_time_in_seconds = 0.0;
357   summary->line_search_total_time_in_seconds  = 0.0;
358   summary->inner_iterations_given             = options.use_inner_iterations;
359   summary->line_search_direction_type         = options.line_search_direction_type;         //  NOLINT
360   summary->line_search_interpolation_type     = options.line_search_interpolation_type;     //  NOLINT
361   summary->line_search_type                   = options.line_search_type;
362   summary->linear_solver_type_given           = options.linear_solver_type;
363   summary->max_lbfgs_rank                     = options.max_lbfgs_rank;
364   summary->minimizer_type                     = options.minimizer_type;
365   summary->nonlinear_conjugate_gradient_type  = options.nonlinear_conjugate_gradient_type;  //  NOLINT
366   summary->num_linear_solver_threads_given    = options.num_linear_solver_threads;          //  NOLINT
367   summary->num_threads_given                  = options.num_threads;
368   summary->preconditioner_type_given          = options.preconditioner_type;
369   summary->sparse_linear_algebra_library_type = options.sparse_linear_algebra_library_type; //  NOLINT
370   summary->trust_region_strategy_type         = options.trust_region_strategy_type;         //  NOLINT
371   summary->visibility_clustering_type         = options.visibility_clustering_type;         //  NOLINT
372 }
373
374 void PostSolveSummarize(const internal::PreprocessedProblem& pp,
375                         Solver::Summary* summary) {
376   internal::OrderingToGroupSizes(pp.options.linear_solver_ordering.get(),
377                                  &(summary->linear_solver_ordering_used));
378   internal::OrderingToGroupSizes(pp.options.inner_iteration_ordering.get(),
379                                  &(summary->inner_iteration_ordering_used));
380
381   summary->inner_iterations_used          = pp.inner_iteration_minimizer.get() != NULL;     // NOLINT
382   summary->linear_solver_type_used        = pp.options.linear_solver_type;
383   summary->num_linear_solver_threads_used = pp.options.num_linear_solver_threads;           // NOLINT
384   summary->num_threads_used               = pp.options.num_threads;
385   summary->preconditioner_type_used       = pp.options.preconditioner_type;                 // NOLINT
386
387   internal::SetSummaryFinalCost(summary);
388
389   if (pp.reduced_program.get() != NULL) {
390     SummarizeReducedProgram(*pp.reduced_program, summary);
391   }
392
393   // It is possible that no evaluator was created. This would be the
394   // case if the preprocessor failed, or if the reduced problem did
395   // not contain any parameter blocks. Thus, only extract the
396   // evaluator statistics if one exists.
397   if (pp.evaluator.get() != NULL) {
398     const map<string, double>& evaluator_time_statistics =
399         pp.evaluator->TimeStatistics();
400     summary->residual_evaluation_time_in_seconds =
401         FindWithDefault(evaluator_time_statistics, "Evaluator::Residual", 0.0);
402     summary->jacobian_evaluation_time_in_seconds =
403         FindWithDefault(evaluator_time_statistics, "Evaluator::Jacobian", 0.0);
404   }
405
406   // Again, like the evaluator, there may or may not be a linear
407   // solver from which we can extract run time statistics. In
408   // particular the line search solver does not use a linear solver.
409   if (pp.linear_solver.get() != NULL) {
410     const map<string, double>& linear_solver_time_statistics =
411         pp.linear_solver->TimeStatistics();
412     summary->linear_solver_time_in_seconds =
413         FindWithDefault(linear_solver_time_statistics,
414                         "LinearSolver::Solve",
415                         0.0);
416   }
417 }
418
419 void Minimize(internal::PreprocessedProblem* pp,
420               Solver::Summary* summary) {
421   using internal::Program;
422   using internal::scoped_ptr;
423   using internal::Minimizer;
424
425   Program* program = pp->reduced_program.get();
426   if (pp->reduced_program->NumParameterBlocks() == 0) {
427     summary->message = "Function tolerance reached. "
428         "No non-constant parameter blocks found.";
429     summary->termination_type = CONVERGENCE;
430     VLOG_IF(1, pp->options.logging_type != SILENT) << summary->message;
431     summary->initial_cost = summary->fixed_cost;
432     summary->final_cost = summary->fixed_cost;
433     return;
434   }
435
436   scoped_ptr<Minimizer> minimizer(
437       Minimizer::Create(pp->options.minimizer_type));
438   minimizer->Minimize(pp->minimizer_options,
439                       pp->reduced_parameters.data(),
440                       summary);
441
442   if (summary->IsSolutionUsable()) {
443     program->StateVectorToParameterBlocks(pp->reduced_parameters.data());
444     program->CopyParameterBlockStateToUserState();
445   }
446 }
447
448 }  // namespace
449
450 bool Solver::Options::IsValid(string* error) const {
451   if (!CommonOptionsAreValid(*this, error)) {
452     return false;
453   }
454
455   if (minimizer_type == TRUST_REGION &&
456       !TrustRegionOptionsAreValid(*this, error)) {
457     return false;
458   }
459
460   // We do not know if the problem is bounds constrained or not, if it
461   // is then the trust region solver will also use the line search
462   // solver to do a projection onto the box constraints, so make sure
463   // that the line search options are checked independent of what
464   // minimizer algorithm is being used.
465   return LineSearchOptionsAreValid(*this, error);
466 }
467
468 Solver::~Solver() {}
469
470 void Solver::Solve(const Solver::Options& options,
471                    Problem* problem,
472                    Solver::Summary* summary) {
473   using internal::PreprocessedProblem;
474   using internal::Preprocessor;
475   using internal::ProblemImpl;
476   using internal::Program;
477   using internal::scoped_ptr;
478   using internal::WallTimeInSeconds;
479
480   CHECK_NOTNULL(problem);
481   CHECK_NOTNULL(summary);
482
483   double start_time = WallTimeInSeconds();
484   *summary = Summary();
485   if (!options.IsValid(&summary->message)) {
486     LOG(ERROR) << "Terminating: " << summary->message;
487     return;
488   }
489
490   ProblemImpl* problem_impl = problem->problem_impl_.get();
491   Program* program = problem_impl->mutable_program();
492   PreSolveSummarize(options, problem_impl, summary);
493
494   // Make sure that all the parameter blocks states are set to the
495   // values provided by the user.
496   program->SetParameterBlockStatePtrsToUserStatePtrs();
497
498   scoped_ptr<internal::ProblemImpl> gradient_checking_problem;
499   if (options.check_gradients) {
500     gradient_checking_problem.reset(
501         CreateGradientCheckingProblemImpl(
502             problem_impl,
503             options.numeric_derivative_relative_step_size,
504             options.gradient_check_relative_precision));
505     problem_impl = gradient_checking_problem.get();
506     program = problem_impl->mutable_program();
507   }
508
509   scoped_ptr<Preprocessor> preprocessor(
510       Preprocessor::Create(options.minimizer_type));
511   PreprocessedProblem pp;
512   const bool status = preprocessor->Preprocess(options, problem_impl, &pp);
513   summary->fixed_cost = pp.fixed_cost;
514   summary->preprocessor_time_in_seconds = WallTimeInSeconds() - start_time;
515
516   if (status) {
517     const double minimizer_start_time = WallTimeInSeconds();
518     Minimize(&pp, summary);
519     summary->minimizer_time_in_seconds =
520         WallTimeInSeconds() - minimizer_start_time;
521   } else {
522     summary->message = pp.error;
523   }
524
525   const double postprocessor_start_time = WallTimeInSeconds();
526   problem_impl = problem->problem_impl_.get();
527   program = problem_impl->mutable_program();
528   // On exit, ensure that the parameter blocks again point at the user
529   // provided values and the parameter blocks are numbered according
530   // to their position in the original user provided program.
531   program->SetParameterBlockStatePtrsToUserStatePtrs();
532   program->SetParameterOffsetsAndIndex();
533   PostSolveSummarize(pp, summary);
534   summary->postprocessor_time_in_seconds =
535       WallTimeInSeconds() - postprocessor_start_time;
536
537   summary->total_time_in_seconds = WallTimeInSeconds() - start_time;
538 }
539
540 void Solve(const Solver::Options& options,
541            Problem* problem,
542            Solver::Summary* summary) {
543   Solver solver;
544   solver.Solve(options, problem, summary);
545 }
546
547 Solver::Summary::Summary()
548     // Invalid values for most fields, to ensure that we are not
549     // accidentally reporting default values.
550     : minimizer_type(TRUST_REGION),
551       termination_type(FAILURE),
552       message("ceres::Solve was not called."),
553       initial_cost(-1.0),
554       final_cost(-1.0),
555       fixed_cost(-1.0),
556       num_successful_steps(-1),
557       num_unsuccessful_steps(-1),
558       num_inner_iteration_steps(-1),
559       preprocessor_time_in_seconds(-1.0),
560       minimizer_time_in_seconds(-1.0),
561       postprocessor_time_in_seconds(-1.0),
562       total_time_in_seconds(-1.0),
563       linear_solver_time_in_seconds(-1.0),
564       residual_evaluation_time_in_seconds(-1.0),
565       jacobian_evaluation_time_in_seconds(-1.0),
566       inner_iteration_time_in_seconds(-1.0),
567       line_search_cost_evaluation_time_in_seconds(-1.0),
568       line_search_gradient_evaluation_time_in_seconds(-1.0),
569       line_search_polynomial_minimization_time_in_seconds(-1.0),
570       line_search_total_time_in_seconds(-1.0),
571       num_parameter_blocks(-1),
572       num_parameters(-1),
573       num_effective_parameters(-1),
574       num_residual_blocks(-1),
575       num_residuals(-1),
576       num_parameter_blocks_reduced(-1),
577       num_parameters_reduced(-1),
578       num_effective_parameters_reduced(-1),
579       num_residual_blocks_reduced(-1),
580       num_residuals_reduced(-1),
581       is_constrained(false),
582       num_threads_given(-1),
583       num_threads_used(-1),
584       num_linear_solver_threads_given(-1),
585       num_linear_solver_threads_used(-1),
586       linear_solver_type_given(SPARSE_NORMAL_CHOLESKY),
587       linear_solver_type_used(SPARSE_NORMAL_CHOLESKY),
588       inner_iterations_given(false),
589       inner_iterations_used(false),
590       preconditioner_type_given(IDENTITY),
591       preconditioner_type_used(IDENTITY),
592       visibility_clustering_type(CANONICAL_VIEWS),
593       trust_region_strategy_type(LEVENBERG_MARQUARDT),
594       dense_linear_algebra_library_type(EIGEN),
595       sparse_linear_algebra_library_type(SUITE_SPARSE),
596       line_search_direction_type(LBFGS),
597       line_search_type(ARMIJO),
598       line_search_interpolation_type(BISECTION),
599       nonlinear_conjugate_gradient_type(FLETCHER_REEVES),
600       max_lbfgs_rank(-1) {
601 }
602
603 using internal::StringAppendF;
604 using internal::StringPrintf;
605
606 string Solver::Summary::BriefReport() const {
607   return StringPrintf("Ceres Solver Report: "
608                       "Iterations: %d, "
609                       "Initial cost: %e, "
610                       "Final cost: %e, "
611                       "Termination: %s",
612                       num_successful_steps + num_unsuccessful_steps,
613                       initial_cost,
614                       final_cost,
615                       TerminationTypeToString(termination_type));
616 }
617
618 string Solver::Summary::FullReport() const {
619   using internal::VersionString;
620
621   string report = string("\nSolver Summary (v " + VersionString() + ")\n\n");
622
623   StringAppendF(&report, "%45s    %21s\n", "Original", "Reduced");
624   StringAppendF(&report, "Parameter blocks    % 25d% 25d\n",
625                 num_parameter_blocks, num_parameter_blocks_reduced);
626   StringAppendF(&report, "Parameters          % 25d% 25d\n",
627                 num_parameters, num_parameters_reduced);
628   if (num_effective_parameters_reduced != num_parameters_reduced) {
629     StringAppendF(&report, "Effective parameters% 25d% 25d\n",
630                   num_effective_parameters, num_effective_parameters_reduced);
631   }
632   StringAppendF(&report, "Residual blocks     % 25d% 25d\n",
633                 num_residual_blocks, num_residual_blocks_reduced);
634   StringAppendF(&report, "Residual            % 25d% 25d\n",
635                 num_residuals, num_residuals_reduced);
636
637   if (minimizer_type == TRUST_REGION) {
638     // TRUST_SEARCH HEADER
639     StringAppendF(&report, "\nMinimizer                 %19s\n",
640                   "TRUST_REGION");
641
642     if (linear_solver_type_used == DENSE_NORMAL_CHOLESKY ||
643         linear_solver_type_used == DENSE_SCHUR ||
644         linear_solver_type_used == DENSE_QR) {
645       StringAppendF(&report, "\nDense linear algebra library  %15s\n",
646                     DenseLinearAlgebraLibraryTypeToString(
647                         dense_linear_algebra_library_type));
648     }
649
650     if (linear_solver_type_used == SPARSE_NORMAL_CHOLESKY ||
651         linear_solver_type_used == SPARSE_SCHUR ||
652         (linear_solver_type_used == ITERATIVE_SCHUR &&
653          (preconditioner_type_used == CLUSTER_JACOBI ||
654           preconditioner_type_used == CLUSTER_TRIDIAGONAL))) {
655       StringAppendF(&report, "\nSparse linear algebra library %15s\n",
656                     SparseLinearAlgebraLibraryTypeToString(
657                         sparse_linear_algebra_library_type));
658     }
659
660     StringAppendF(&report, "Trust region strategy     %19s",
661                   TrustRegionStrategyTypeToString(
662                       trust_region_strategy_type));
663     if (trust_region_strategy_type == DOGLEG) {
664       if (dogleg_type == TRADITIONAL_DOGLEG) {
665         StringAppendF(&report, " (TRADITIONAL)");
666       } else {
667         StringAppendF(&report, " (SUBSPACE)");
668       }
669     }
670     StringAppendF(&report, "\n");
671     StringAppendF(&report, "\n");
672
673     StringAppendF(&report, "%45s    %21s\n", "Given",  "Used");
674     StringAppendF(&report, "Linear solver       %25s%25s\n",
675                   LinearSolverTypeToString(linear_solver_type_given),
676                   LinearSolverTypeToString(linear_solver_type_used));
677
678     if (linear_solver_type_given == CGNR ||
679         linear_solver_type_given == ITERATIVE_SCHUR) {
680       StringAppendF(&report, "Preconditioner      %25s%25s\n",
681                     PreconditionerTypeToString(preconditioner_type_given),
682                     PreconditionerTypeToString(preconditioner_type_used));
683     }
684
685     if (preconditioner_type_used == CLUSTER_JACOBI ||
686         preconditioner_type_used == CLUSTER_TRIDIAGONAL) {
687       StringAppendF(&report, "Visibility clustering%24s%25s\n",
688                     VisibilityClusteringTypeToString(
689                         visibility_clustering_type),
690                     VisibilityClusteringTypeToString(
691                         visibility_clustering_type));
692     }
693     StringAppendF(&report, "Threads             % 25d% 25d\n",
694                   num_threads_given, num_threads_used);
695     StringAppendF(&report, "Linear solver threads % 23d% 25d\n",
696                   num_linear_solver_threads_given,
697                   num_linear_solver_threads_used);
698
699     if (IsSchurType(linear_solver_type_used)) {
700       string given;
701       StringifyOrdering(linear_solver_ordering_given, &given);
702       string used;
703       StringifyOrdering(linear_solver_ordering_used, &used);
704       StringAppendF(&report,
705                     "Linear solver ordering %22s %24s\n",
706                     given.c_str(),
707                     used.c_str());
708     }
709
710     if (inner_iterations_given) {
711       StringAppendF(&report,
712                     "Use inner iterations     %20s     %20s\n",
713                     inner_iterations_given ? "True" : "False",
714                     inner_iterations_used ? "True" : "False");
715     }
716
717     if (inner_iterations_used) {
718       string given;
719       StringifyOrdering(inner_iteration_ordering_given, &given);
720       string used;
721       StringifyOrdering(inner_iteration_ordering_used, &used);
722     StringAppendF(&report,
723                   "Inner iteration ordering %20s %24s\n",
724                   given.c_str(),
725                   used.c_str());
726     }
727   } else {
728     // LINE_SEARCH HEADER
729     StringAppendF(&report, "\nMinimizer                 %19s\n", "LINE_SEARCH");
730
731
732     string line_search_direction_string;
733     if (line_search_direction_type == LBFGS) {
734       line_search_direction_string = StringPrintf("LBFGS (%d)", max_lbfgs_rank);
735     } else if (line_search_direction_type == NONLINEAR_CONJUGATE_GRADIENT) {
736       line_search_direction_string =
737           NonlinearConjugateGradientTypeToString(
738               nonlinear_conjugate_gradient_type);
739     } else {
740       line_search_direction_string =
741           LineSearchDirectionTypeToString(line_search_direction_type);
742     }
743
744     StringAppendF(&report, "Line search direction     %19s\n",
745                   line_search_direction_string.c_str());
746
747     const string line_search_type_string =
748         StringPrintf("%s %s",
749                      LineSearchInterpolationTypeToString(
750                          line_search_interpolation_type),
751                      LineSearchTypeToString(line_search_type));
752     StringAppendF(&report, "Line search type          %19s\n",
753                   line_search_type_string.c_str());
754     StringAppendF(&report, "\n");
755
756     StringAppendF(&report, "%45s    %21s\n", "Given",  "Used");
757     StringAppendF(&report, "Threads             % 25d% 25d\n",
758                   num_threads_given, num_threads_used);
759   }
760
761   StringAppendF(&report, "\nCost:\n");
762   StringAppendF(&report, "Initial        % 30e\n", initial_cost);
763   if (termination_type != FAILURE &&
764       termination_type != USER_FAILURE) {
765     StringAppendF(&report, "Final          % 30e\n", final_cost);
766     StringAppendF(&report, "Change         % 30e\n",
767                   initial_cost - final_cost);
768   }
769
770   StringAppendF(&report, "\nMinimizer iterations         % 16d\n",
771                 num_successful_steps + num_unsuccessful_steps);
772
773   // Successful/Unsuccessful steps only matter in the case of the
774   // trust region solver. Line search terminates when it encounters
775   // the first unsuccessful step.
776   if (minimizer_type == TRUST_REGION) {
777     StringAppendF(&report, "Successful steps               % 14d\n",
778                   num_successful_steps);
779     StringAppendF(&report, "Unsuccessful steps             % 14d\n",
780                   num_unsuccessful_steps);
781   }
782   if (inner_iterations_used) {
783     StringAppendF(&report, "Steps with inner iterations    % 14d\n",
784                   num_inner_iteration_steps);
785   }
786
787   const bool print_line_search_timing_information =
788       minimizer_type == LINE_SEARCH ||
789       (minimizer_type == TRUST_REGION && is_constrained);
790
791   StringAppendF(&report, "\nTime (in seconds):\n");
792   StringAppendF(&report, "Preprocessor        %25.4f\n",
793                 preprocessor_time_in_seconds);
794
795   StringAppendF(&report, "\n  Residual evaluation %23.4f\n",
796                 residual_evaluation_time_in_seconds);
797   if (print_line_search_timing_information) {
798     StringAppendF(&report, "    Line search cost evaluation    %10.4f\n",
799                   line_search_cost_evaluation_time_in_seconds);
800   }
801   StringAppendF(&report, "  Jacobian evaluation %23.4f\n",
802                 jacobian_evaluation_time_in_seconds);
803   if (print_line_search_timing_information) {
804     StringAppendF(&report, "    Line search gradient evaluation    %6.4f\n",
805                   line_search_gradient_evaluation_time_in_seconds);
806   }
807
808   if (minimizer_type == TRUST_REGION) {
809     StringAppendF(&report, "  Linear solver       %23.4f\n",
810                   linear_solver_time_in_seconds);
811   }
812
813   if (inner_iterations_used) {
814     StringAppendF(&report, "  Inner iterations    %23.4f\n",
815                   inner_iteration_time_in_seconds);
816   }
817
818   if (print_line_search_timing_information) {
819     StringAppendF(&report, "  Line search polynomial minimization  %.4f\n",
820                   line_search_polynomial_minimization_time_in_seconds);
821   }
822
823   StringAppendF(&report, "Minimizer           %25.4f\n\n",
824                 minimizer_time_in_seconds);
825
826   StringAppendF(&report, "Postprocessor        %24.4f\n",
827                 postprocessor_time_in_seconds);
828
829   StringAppendF(&report, "Total               %25.4f\n\n",
830                 total_time_in_seconds);
831
832   StringAppendF(&report, "Termination:        %25s (%s)\n",
833                 TerminationTypeToString(termination_type), message.c_str());
834   return report;
835 }
836
837 bool Solver::Summary::IsSolutionUsable() const {
838   return internal::IsSolutionUsable(*this);
839 }
840
841 }  // namespace ceres