internal/ceres/solver_test.cc

   1 // Ceres Solver - A fast non-linear least squares minimizer
   2 // Copyright 2015 Google Inc. All rights reserved.
   3 // http://ceres-solver.org/
   4 //
   5 // Redistribution and use in source and binary forms, with or without
   6 // modification, are permitted provided that the following conditions are met:
   7 //
   8 // * Redistributions of source code must retain the above copyright notice,
   9 //   this list of conditions and the following disclaimer.
  10 // * Redistributions in binary form must reproduce the above copyright notice,
  11 //   this list of conditions and the following disclaimer in the documentation
  12 //   and/or other materials provided with the distribution.
  13 // * Neither the name of Google Inc. nor the names of its contributors may be
  14 //   used to endorse or promote products derived from this software without
  15 //   specific prior written permission.
  16 //
  17 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  18 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  19 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  20 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  21 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  22 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  23 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  24 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  25 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  26 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  27 // POSSIBILITY OF SUCH DAMAGE.
  28 //
  29 // Author: sameeragarwal@google.com (Sameer Agarwal)
  30
  31 #include "ceres/solver.h"
  32
  33 #include <limits>
  34 #include <cmath>
  35 #include <vector>
  36 #include "gtest/gtest.h"
  37 #include "ceres/internal/scoped_ptr.h"
  38 #include "ceres/autodiff_cost_function.h"
  39 #include "ceres/sized_cost_function.h"
  40 #include "ceres/problem.h"
  41 #include "ceres/problem_impl.h"
  42
  43 namespace ceres {
  44 namespace internal {
  45
  46 using std::string;
  47
  48 TEST(SolverOptions, DefaultTrustRegionOptionsAreValid) {
  49   Solver::Options options;
  50   options.minimizer_type = TRUST_REGION;
  51   string error;
  52   EXPECT_TRUE(options.IsValid(&error)) << error;
  53 }
  54
  55 TEST(SolverOptions, DefaultLineSearchOptionsAreValid) {
  56   Solver::Options options;
  57   options.minimizer_type = LINE_SEARCH;
  58   string error;
  59   EXPECT_TRUE(options.IsValid(&error)) << error;
  60 }
  61
  62 struct QuadraticCostFunctor {
  63   template <typename T> bool operator()(const T* const x,
  64                                         T* residual) const {
  65     residual[0] = T(5.0) - *x;
  66     return true;
  67   }
  68
  69   static CostFunction* Create() {
  70     return new AutoDiffCostFunction<QuadraticCostFunctor, 1, 1>(
  71         new QuadraticCostFunctor);
  72   }
  73 };
  74
  75 struct RememberingCallback : public IterationCallback {
  76   explicit RememberingCallback(double *x) : calls(0), x(x) {}
  77   virtual ~RememberingCallback() {}
  78   virtual CallbackReturnType operator()(const IterationSummary& summary) {
  79     x_values.push_back(*x);
  80     return SOLVER_CONTINUE;
  81   }
  82   int calls;
  83   double *x;
  84   std::vector<double> x_values;
  85 };
  86
  87 TEST(Solver, UpdateStateEveryIterationOption) {
  88   double x = 50.0;
  89   const double original_x = x;
  90
  91   scoped_ptr<CostFunction> cost_function(QuadraticCostFunctor::Create());
  92   Problem::Options problem_options;
  93   problem_options.cost_function_ownership = DO_NOT_TAKE_OWNERSHIP;
  94   Problem problem(problem_options);
  95   problem.AddResidualBlock(cost_function.get(), NULL, &x);
  96
  97   Solver::Options options;
  98   options.linear_solver_type = DENSE_QR;
  99
 100   RememberingCallback callback(&x);
 101   options.callbacks.push_back(&callback);
 102
 103   Solver::Summary summary;
 104
 105   int num_iterations;
 106
 107   // First try: no updating.
 108   Solve(options, &problem, &summary);
 109   num_iterations = summary.num_successful_steps +
 110                    summary.num_unsuccessful_steps;
 111   EXPECT_GT(num_iterations, 1);
 112   for (int i = 0; i < callback.x_values.size(); ++i) {
 113     EXPECT_EQ(50.0, callback.x_values[i]);
 114   }
 115
 116   // Second try: with updating
 117   x = 50.0;
 118   options.update_state_every_iteration = true;
 119   callback.x_values.clear();
 120   Solve(options, &problem, &summary);
 121   num_iterations = summary.num_successful_steps +
 122                    summary.num_unsuccessful_steps;
 123   EXPECT_GT(num_iterations, 1);
 124   EXPECT_EQ(original_x, callback.x_values[0]);
 125   EXPECT_NE(original_x, callback.x_values[1]);
 126 }
 127
 128 // The parameters must be in separate blocks so that they can be individually
 129 // set constant or not.
 130 struct Quadratic4DCostFunction {
 131   template <typename T> bool operator()(const T* const x,
 132                                         const T* const y,
 133                                         const T* const z,
 134                                         const T* const w,
 135                                         T* residual) const {
 136     // A 4-dimension axis-aligned quadratic.
 137     residual[0] = T(10.0) - *x +
 138                   T(20.0) - *y +
 139                   T(30.0) - *z +
 140                   T(40.0) - *w;
 141     return true;
 142   }
 143
 144   static CostFunction* Create() {
 145     return new AutoDiffCostFunction<Quadratic4DCostFunction, 1, 1, 1, 1, 1>(
 146         new Quadratic4DCostFunction);
 147   }
 148 };
 149
 150 // A cost function that simply returns its argument.
 151 class UnaryIdentityCostFunction : public SizedCostFunction<1, 1> {
 152  public:
 153   virtual bool Evaluate(double const* const* parameters,
 154                         double* residuals,
 155                         double** jacobians) const {
 156     residuals[0] = parameters[0][0];
 157     if (jacobians != NULL && jacobians[0] != NULL) {
 158       jacobians[0][0] = 1.0;
 159     }
 160     return true;
 161   }
 162 };
 163
 164 TEST(Solver, TrustRegionProblemHasNoParameterBlocks) {
 165   Problem problem;
 166   Solver::Options options;
 167   options.minimizer_type = TRUST_REGION;
 168   Solver::Summary summary;
 169   Solve(options, &problem, &summary);
 170   EXPECT_EQ(summary.termination_type, CONVERGENCE);
 171   EXPECT_EQ(summary.message,
 172             "Function tolerance reached. "
 173             "No non-constant parameter blocks found.");
 174 }
 175
 176 TEST(Solver, LineSearchProblemHasNoParameterBlocks) {
 177   Problem problem;
 178   Solver::Options options;
 179   options.minimizer_type = LINE_SEARCH;
 180   Solver::Summary summary;
 181   Solve(options, &problem, &summary);
 182   EXPECT_EQ(summary.termination_type, CONVERGENCE);
 183   EXPECT_EQ(summary.message,
 184             "Function tolerance reached. "
 185             "No non-constant parameter blocks found.");
 186 }
 187
 188 TEST(Solver, TrustRegionProblemHasZeroResiduals) {
 189   Problem problem;
 190   double x = 1;
 191   problem.AddParameterBlock(&x, 1);
 192   Solver::Options options;
 193   options.minimizer_type = TRUST_REGION;
 194   Solver::Summary summary;
 195   Solve(options, &problem, &summary);
 196   EXPECT_EQ(summary.termination_type, CONVERGENCE);
 197   EXPECT_EQ(summary.message,
 198             "Function tolerance reached. "
 199             "No non-constant parameter blocks found.");
 200 }
 201
 202 TEST(Solver, LineSearchProblemHasZeroResiduals) {
 203   Problem problem;
 204   double x = 1;
 205   problem.AddParameterBlock(&x, 1);
 206   Solver::Options options;
 207   options.minimizer_type = LINE_SEARCH;
 208   Solver::Summary summary;
 209   Solve(options, &problem, &summary);
 210   EXPECT_EQ(summary.termination_type, CONVERGENCE);
 211   EXPECT_EQ(summary.message,
 212             "Function tolerance reached. "
 213             "No non-constant parameter blocks found.");
 214 }
 215
 216 TEST(Solver, TrustRegionProblemIsConstant) {
 217   Problem problem;
 218   double x = 1;
 219   problem.AddResidualBlock(new UnaryIdentityCostFunction, NULL, &x);
 220   problem.SetParameterBlockConstant(&x);
 221   Solver::Options options;
 222   options.minimizer_type = TRUST_REGION;
 223   Solver::Summary summary;
 224   Solve(options, &problem, &summary);
 225   EXPECT_EQ(summary.termination_type, CONVERGENCE);
 226   EXPECT_EQ(summary.initial_cost, 1.0 / 2.0);
 227   EXPECT_EQ(summary.final_cost, 1.0 / 2.0);
 228 }
 229
 230 TEST(Solver, LineSearchProblemIsConstant) {
 231   Problem problem;
 232   double x = 1;
 233   problem.AddResidualBlock(new UnaryIdentityCostFunction, NULL, &x);
 234   problem.SetParameterBlockConstant(&x);
 235   Solver::Options options;
 236   options.minimizer_type = LINE_SEARCH;
 237   Solver::Summary summary;
 238   Solve(options, &problem, &summary);
 239   EXPECT_EQ(summary.termination_type, CONVERGENCE);
 240   EXPECT_EQ(summary.initial_cost, 1.0 / 2.0);
 241   EXPECT_EQ(summary.final_cost, 1.0 / 2.0);
 242 }
 243
 244 #if defined(CERES_NO_SUITESPARSE)
 245 TEST(Solver, SparseNormalCholeskyNoSuiteSparse) {
 246   Solver::Options options;
 247   options.sparse_linear_algebra_library_type = SUITE_SPARSE;
 248   options.linear_solver_type = SPARSE_NORMAL_CHOLESKY;
 249   string message;
 250   EXPECT_FALSE(options.IsValid(&message));
 251 }
 252
 253 TEST(Solver, SparseSchurNoSuiteSparse) {
 254   Solver::Options options;
 255   options.sparse_linear_algebra_library_type = SUITE_SPARSE;
 256   options.linear_solver_type = SPARSE_SCHUR;
 257   string message;
 258   EXPECT_FALSE(options.IsValid(&message));
 259 }
 260 #endif
 261
 262 #if defined(CERES_NO_CXSPARSE)
 263 TEST(Solver, SparseNormalCholeskyNoCXSparse) {
 264   Solver::Options options;
 265   options.sparse_linear_algebra_library_type = CX_SPARSE;
 266   options.linear_solver_type = SPARSE_NORMAL_CHOLESKY;
 267   string message;
 268   EXPECT_FALSE(options.IsValid(&message));
 269 }
 270
 271 TEST(Solver, SparseSchurNoCXSparse) {
 272   Solver::Options options;
 273   options.sparse_linear_algebra_library_type = CX_SPARSE;
 274   options.linear_solver_type = SPARSE_SCHUR;
 275   string message;
 276   EXPECT_FALSE(options.IsValid(&message));
 277 }
 278 #endif
 279
 280 #if !defined(CERES_USE_EIGEN_SPARSE)
 281 TEST(Solver, SparseNormalCholeskyNoEigenSparse) {
 282   Solver::Options options;
 283   options.sparse_linear_algebra_library_type = EIGEN_SPARSE;
 284   options.linear_solver_type = SPARSE_NORMAL_CHOLESKY;
 285   string message;
 286   EXPECT_FALSE(options.IsValid(&message));
 287 }
 288
 289 TEST(Solver, SparseSchurNoEigenSparse) {
 290   Solver::Options options;
 291   options.sparse_linear_algebra_library_type = EIGEN_SPARSE;
 292   options.linear_solver_type = SPARSE_SCHUR;
 293   string message;
 294   EXPECT_FALSE(options.IsValid(&message));
 295 }
 296 #endif
 297
 298 TEST(Solver, SparseNormalCholeskyNoSparseLibrary) {
 299   Solver::Options options;
 300   options.sparse_linear_algebra_library_type = NO_SPARSE;
 301   options.linear_solver_type = SPARSE_NORMAL_CHOLESKY;
 302   string message;
 303   EXPECT_FALSE(options.IsValid(&message));
 304 }
 305
 306 TEST(Solver, SparseSchurNoSparseLibrary) {
 307   Solver::Options options;
 308   options.sparse_linear_algebra_library_type = NO_SPARSE;
 309   options.linear_solver_type = SPARSE_SCHUR;
 310   string message;
 311   EXPECT_FALSE(options.IsValid(&message));
 312 }
 313
 314 TEST(Solver, IterativeSchurWithClusterJacobiPerconditionerNoSparseLibrary) {
 315   Solver::Options options;
 316   options.sparse_linear_algebra_library_type = NO_SPARSE;
 317   options.linear_solver_type = ITERATIVE_SCHUR;
 318   // Requires SuiteSparse.
 319   options.preconditioner_type = CLUSTER_JACOBI;
 320   string message;
 321   EXPECT_FALSE(options.IsValid(&message));
 322 }
 323
 324 TEST(Solver, IterativeSchurWithClusterTridiagonalPerconditionerNoSparseLibrary) {
 325   Solver::Options options;
 326   options.sparse_linear_algebra_library_type = NO_SPARSE;
 327   options.linear_solver_type = ITERATIVE_SCHUR;
 328   // Requires SuiteSparse.
 329   options.preconditioner_type = CLUSTER_TRIDIAGONAL;
 330   string message;
 331   EXPECT_FALSE(options.IsValid(&message));
 332 }
 333
 334 TEST(Solver, IterativeLinearSolverForDogleg) {
 335   Solver::Options options;
 336   options.trust_region_strategy_type = DOGLEG;
 337   string message;
 338   options.linear_solver_type = ITERATIVE_SCHUR;
 339   EXPECT_FALSE(options.IsValid(&message));
 340
 341   options.linear_solver_type = CGNR;
 342   EXPECT_FALSE(options.IsValid(&message));
 343 }
 344
 345 TEST(Solver, LinearSolverTypeNormalOperation) {
 346   Solver::Options options;
 347   options.linear_solver_type = DENSE_QR;
 348
 349   string message;
 350   EXPECT_TRUE(options.IsValid(&message));
 351
 352   options.linear_solver_type = DENSE_NORMAL_CHOLESKY;
 353   EXPECT_TRUE(options.IsValid(&message));
 354
 355   options.linear_solver_type = DENSE_SCHUR;
 356   EXPECT_TRUE(options.IsValid(&message));
 357
 358   options.linear_solver_type = SPARSE_SCHUR;
 359 #if defined(CERES_NO_SUITESPARSE) &&            \
 360     defined(CERES_NO_CXSPARSE) &&               \
 361    !defined(CERES_USE_EIGEN_SPARSE)
 362   EXPECT_FALSE(options.IsValid(&message));
 363 #else
 364   EXPECT_TRUE(options.IsValid(&message));
 365 #endif
 366
 367   options.linear_solver_type = ITERATIVE_SCHUR;
 368   EXPECT_TRUE(options.IsValid(&message));
 369 }
 370
 371 template<int kNumResiduals, int N1 = 0, int N2 = 0, int N3 = 0>
 372 class DummyCostFunction : public SizedCostFunction<kNumResiduals, N1, N2, N3> {
 373  public:
 374   bool Evaluate(double const* const* parameters,
 375                 double* residuals,
 376                 double** jacobians) const {
 377     for (int i = 0; i < kNumResiduals; ++i) {
 378       residuals[i] = kNumResiduals * kNumResiduals + i;
 379     }
 380
 381     return true;
 382   }
 383 };
 384
 385 TEST(Solver, FixedCostForConstantProblem) {
 386   double x = 1.0;
 387   Problem problem;
 388   problem.AddResidualBlock(new DummyCostFunction<2, 1>(), NULL, &x);
 389   problem.SetParameterBlockConstant(&x);
 390   const double expected_cost = 41.0 / 2.0;  // 1/2 * ((4 + 0)^2 + (4 + 1)^2)
 391   Solver::Options options;
 392   Solver::Summary summary;
 393   Solve(options, &problem, &summary);
 394   EXPECT_TRUE(summary.IsSolutionUsable());
 395   EXPECT_EQ(summary.fixed_cost, expected_cost);
 396   EXPECT_EQ(summary.initial_cost, expected_cost);
 397   EXPECT_EQ(summary.final_cost, expected_cost);
 398   EXPECT_EQ(summary.iterations.size(), 0);
 399 }
 400
 401 }  // namespace internal
 402 }  // namespace ceres