internal/ceres/numeric_diff_test_utils.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 //         tbennun@gmail.com (Tal Ben-Nun)
  31
  32 #include "ceres/numeric_diff_test_utils.h"
  33
  34 #include <algorithm>
  35 #include <cmath>
  36 #include "ceres/cost_function.h"
  37 #include "ceres/internal/macros.h"
  38 #include "ceres/test_util.h"
  39 #include "ceres/types.h"
  40 #include "gtest/gtest.h"
  41
  42
  43 namespace ceres {
  44 namespace internal {
  45
  46 bool EasyFunctor::operator()(const double* x1,
  47                              const double* x2,
  48                              double* residuals) const {
  49   residuals[0] = residuals[1] = residuals[2] = 0;
  50   for (int i = 0; i < 5; ++i) {
  51     residuals[0] += x1[i] * x2[i];
  52     residuals[2] += x2[i] * x2[i];
  53   }
  54   residuals[1] = residuals[0] * residuals[0];
  55   return true;
  56 }
  57
  58 void EasyFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect(
  59     const CostFunction& cost_function,
  60     NumericDiffMethodType method) const {
  61   // The x1[0] is made deliberately small to test the performance near
  62   // zero.
  63   double x1[] = { 1e-64, 2.0, 3.0, 4.0, 5.0 };
  64   double x2[] = { 9.0, 9.0, 5.0, 5.0, 1.0 };
  65   double *parameters[] = { &x1[0], &x2[0] };
  66
  67   double dydx1[15];  // 3 x 5, row major.
  68   double dydx2[15];  // 3 x 5, row major.
  69   double *jacobians[2] = { &dydx1[0], &dydx2[0] };
  70
  71   double residuals[3] = {-1e-100, -2e-100, -3e-100 };
  72
  73   ASSERT_TRUE(cost_function.Evaluate(&parameters[0],
  74                                      &residuals[0],
  75                                      &jacobians[0]));
  76
  77   double expected_residuals[3];
  78   EasyFunctor functor;
  79   functor(x1, x2, expected_residuals);
  80   EXPECT_EQ(expected_residuals[0], residuals[0]);
  81   EXPECT_EQ(expected_residuals[1], residuals[1]);
  82   EXPECT_EQ(expected_residuals[2], residuals[2]);
  83
  84   double tolerance = 0.0;
  85   switch (method) {
  86     default:
  87     case CENTRAL:
  88       tolerance = 3e-9;
  89       break;
  90
  91     case FORWARD:
  92       tolerance = 2e-5;
  93       break;
  94
  95     case RIDDERS:
  96       tolerance = 1e-13;
  97       break;
  98   }
  99
 100   for (int i = 0; i < 5; ++i) {
 101     ExpectClose(x2[i],                    dydx1[5 * 0 + i], tolerance);  // y1
 102     ExpectClose(x1[i],                    dydx2[5 * 0 + i], tolerance);
 103     ExpectClose(2 * x2[i] * residuals[0], dydx1[5 * 1 + i], tolerance);  // y2
 104     ExpectClose(2 * x1[i] * residuals[0], dydx2[5 * 1 + i], tolerance);
 105     ExpectClose(0.0,                      dydx1[5 * 2 + i], tolerance);  // y3
 106     ExpectClose(2 * x2[i],                dydx2[5 * 2 + i], tolerance);
 107   }
 108 }
 109
 110 bool TranscendentalFunctor::operator()(const double* x1,
 111                                        const double* x2,
 112                                        double* residuals) const {
 113   double x1x2 = 0;
 114   for (int i = 0; i < 5; ++i) {
 115     x1x2 += x1[i] * x2[i];
 116   }
 117   residuals[0] = sin(x1x2);
 118   residuals[1] = exp(-x1x2 / 10);
 119   return true;
 120 }
 121
 122 void TranscendentalFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect(
 123     const CostFunction& cost_function,
 124     NumericDiffMethodType method) const {
 125   struct {
 126     double x1[5];
 127     double x2[5];
 128   } kTests[] = {
 129     { { 1.0, 2.0, 3.0, 4.0, 5.0 },  // No zeros.
 130       { 9.0, 9.0, 5.0, 5.0, 1.0 },
 131     },
 132     { { 0.0, 2.0, 3.0, 0.0, 5.0 },  // Some zeros x1.
 133       { 9.0, 9.0, 5.0, 5.0, 1.0 },
 134     },
 135     { { 1.0, 2.0, 3.0, 1.0, 5.0 },  // Some zeros x2.
 136       { 0.0, 9.0, 0.0, 5.0, 0.0 },
 137     },
 138     { { 0.0, 0.0, 0.0, 0.0, 0.0 },  // All zeros x1.
 139       { 9.0, 9.0, 5.0, 5.0, 1.0 },
 140     },
 141     { { 1.0, 2.0, 3.0, 4.0, 5.0 },  // All zeros x2.
 142       { 0.0, 0.0, 0.0, 0.0, 0.0 },
 143     },
 144     { { 0.0, 0.0, 0.0, 0.0, 0.0 },  // All zeros.
 145       { 0.0, 0.0, 0.0, 0.0, 0.0 },
 146     },
 147   };
 148
 149   for (int k = 0; k < CERES_ARRAYSIZE(kTests); ++k) {
 150     double *x1 = &(kTests[k].x1[0]);
 151     double *x2 = &(kTests[k].x2[0]);
 152     double *parameters[] = { x1, x2 };
 153
 154     double dydx1[10];
 155     double dydx2[10];
 156     double *jacobians[2] = { &dydx1[0], &dydx2[0] };
 157
 158     double residuals[2];
 159
 160     ASSERT_TRUE(cost_function.Evaluate(&parameters[0],
 161                                        &residuals[0],
 162                                        &jacobians[0]));
 163     double x1x2 = 0;
 164     for (int i = 0; i < 5; ++i) {
 165       x1x2 += x1[i] * x2[i];
 166     }
 167
 168     double tolerance = 0.0;
 169     switch (method) {
 170       default:
 171       case CENTRAL:
 172         tolerance = 2e-7;
 173         break;
 174
 175       case FORWARD:
 176         tolerance = 2e-5;
 177         break;
 178
 179       case RIDDERS:
 180         tolerance = 3e-12;
 181         break;
 182     }
 183
 184     for (int i = 0; i < 5; ++i) {
 185       ExpectClose( x2[i] * cos(x1x2),              dydx1[5 * 0 + i], tolerance);
 186       ExpectClose( x1[i] * cos(x1x2),              dydx2[5 * 0 + i], tolerance);
 187       ExpectClose(-x2[i] * exp(-x1x2 / 10.) / 10., dydx1[5 * 1 + i], tolerance);
 188       ExpectClose(-x1[i] * exp(-x1x2 / 10.) / 10., dydx2[5 * 1 + i], tolerance);
 189     }
 190   }
 191 }
 192
 193 bool ExponentialFunctor::operator()(const double* x1,
 194                                     double* residuals) const {
 195   residuals[0] = exp(x1[0]);
 196   return true;
 197 }
 198
 199 void ExponentialFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect(
 200     const CostFunction& cost_function) const {
 201   // Evaluating the functor at specific points for testing.
 202   double kTests[] = { 1.0, 2.0, 3.0, 4.0, 5.0 };
 203
 204   // Minimal tolerance w.r.t. the cost function and the tests.
 205   const double kTolerance = 2e-14;
 206
 207   for (int k = 0; k < CERES_ARRAYSIZE(kTests); ++k) {
 208     double *parameters[] = { &kTests[k] };
 209     double dydx;
 210     double *jacobians[1] = { &dydx };
 211     double residual;
 212
 213     ASSERT_TRUE(cost_function.Evaluate(&parameters[0],
 214                                        &residual,
 215                                        &jacobians[0]));
 216
 217
 218     double expected_result = exp(kTests[k]);
 219
 220     // Expect residual to be close to exp(x).
 221     ExpectClose(residual, expected_result, kTolerance);
 222
 223     // Check evaluated differences. dydx should also be close to exp(x).
 224     ExpectClose(dydx, expected_result, kTolerance);
 225   }
 226 }
 227
 228 bool RandomizedFunctor::operator()(const double* x1,
 229                                    double* residuals) const {
 230   double random_value = static_cast<double>(rand()) /
 231       static_cast<double>(RAND_MAX);
 232
 233   // Normalize noise to [-factor, factor].
 234   random_value *= 2.0;
 235   random_value -= 1.0;
 236   random_value *= noise_factor_;
 237
 238   residuals[0] = x1[0] * x1[0] + random_value;
 239   return true;
 240 }
 241
 242 void RandomizedFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect(
 243     const CostFunction& cost_function) const {
 244   double kTests[] = { 0.0, 1.0, 3.0, 4.0, 50.0 };
 245
 246   const double kTolerance = 2e-4;
 247
 248   // Initialize random number generator with given seed.
 249   srand(random_seed_);
 250
 251   for (int k = 0; k < CERES_ARRAYSIZE(kTests); ++k) {
 252     double *parameters[] = { &kTests[k] };
 253     double dydx;
 254     double *jacobians[1] = { &dydx };
 255     double residual;
 256
 257     ASSERT_TRUE(cost_function.Evaluate(&parameters[0],
 258                                        &residual,
 259                                        &jacobians[0]));
 260
 261     // Expect residual to be close to x^2 w.r.t. noise factor.
 262     ExpectClose(residual, kTests[k] * kTests[k], noise_factor_);
 263
 264     // Check evaluated differences. (dy/dx = ~2x)
 265     ExpectClose(dydx, 2 * kTests[k], kTolerance);
 266   }
 267 }
 268
 269 }  // namespace internal
 270 }  // namespace ceres