1// Ceres Solver - A fast non-linear least squares minimizer
2// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
3// http://code.google.com/p/ceres-solver/
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: wjr@google.com (William Rucklidge)
30//
31// This file contains tests for the GradientChecker class.
32
33#include "ceres/gradient_checker.h"
34
35#include <cmath>
36#include <cstdlib>
37#include <vector>
38
39#include "ceres/cost_function.h"
40#include "ceres/random.h"
41#include "glog/logging.h"
42#include "gtest/gtest.h"
43
44namespace ceres {
45namespace internal {
46
47// We pick a (non-quadratic) function whose derivative are easy:
48//
49//    f = exp(- a' x).
50//   df = - f a.
51//
52// where 'a' is a vector of the same size as 'x'. In the block
53// version, they are both block vectors, of course.
54class GoodTestTerm : public CostFunction {
55 public:
56  GoodTestTerm(int arity, int const *dim) : arity_(arity) {
57    // Make 'arity' random vectors.
58    a_.resize(arity_);
59    for (int j = 0; j < arity_; ++j) {
60      a_[j].resize(dim[j]);
61      for (int u = 0; u < dim[j]; ++u) {
62        a_[j][u] = 2.0 * RandDouble() - 1.0;
63      }
64    }
65
66    for (int i = 0; i < arity_; i++) {
67      mutable_parameter_block_sizes()->push_back(dim[i]);
68    }
69    set_num_residuals(1);
70  }
71
72  bool Evaluate(double const* const* parameters,
73                double* residuals,
74                double** jacobians) const {
75    // Compute a . x.
76    double ax = 0;
77    for (int j = 0; j < arity_; ++j) {
78      for (int u = 0; u < parameter_block_sizes()[j]; ++u) {
79        ax += a_[j][u] * parameters[j][u];
80      }
81    }
82
83    // This is the cost, but also appears as a factor
84    // in the derivatives.
85    double f = *residuals = exp(-ax);
86
87    // Accumulate 1st order derivatives.
88    if (jacobians) {
89      for (int j = 0; j < arity_; ++j) {
90        if (jacobians[j]) {
91          for (int u = 0; u < parameter_block_sizes()[j]; ++u) {
92            // See comments before class.
93            jacobians[j][u] = - f * a_[j][u];
94          }
95        }
96      }
97    }
98
99    return true;
100  }
101
102 private:
103  int arity_;
104  vector<vector<double> > a_;  // our vectors.
105};
106
107class BadTestTerm : public CostFunction {
108 public:
109  BadTestTerm(int arity, int const *dim) : arity_(arity) {
110    // Make 'arity' random vectors.
111    a_.resize(arity_);
112    for (int j = 0; j < arity_; ++j) {
113      a_[j].resize(dim[j]);
114      for (int u = 0; u < dim[j]; ++u) {
115        a_[j][u] = 2.0 * RandDouble() - 1.0;
116      }
117    }
118
119    for (int i = 0; i < arity_; i++) {
120      mutable_parameter_block_sizes()->push_back(dim[i]);
121    }
122    set_num_residuals(1);
123  }
124
125  bool Evaluate(double const* const* parameters,
126                double* residuals,
127                double** jacobians) const {
128    // Compute a . x.
129    double ax = 0;
130    for (int j = 0; j < arity_; ++j) {
131      for (int u = 0; u < parameter_block_sizes()[j]; ++u) {
132        ax += a_[j][u] * parameters[j][u];
133      }
134    }
135
136    // This is the cost, but also appears as a factor
137    // in the derivatives.
138    double f = *residuals = exp(-ax);
139
140    // Accumulate 1st order derivatives.
141    if (jacobians) {
142      for (int j = 0; j < arity_; ++j) {
143        if (jacobians[j]) {
144          for (int u = 0; u < parameter_block_sizes()[j]; ++u) {
145            // See comments before class.
146            jacobians[j][u] = - f * a_[j][u] + 0.001;
147          }
148        }
149      }
150    }
151
152    return true;
153  }
154
155 private:
156  int arity_;
157  vector<vector<double> > a_;  // our vectors.
158};
159
160TEST(GradientChecker, SmokeTest) {
161  srand(5);
162
163  // Test with 3 blocks of size 2, 3 and 4.
164  int const arity = 3;
165  int const dim[arity] = { 2, 3, 4 };
166
167  // Make a random set of blocks.
168  FixedArray<double*> parameters(arity);
169  for (int j = 0; j < arity; ++j) {
170    parameters[j] = new double[dim[j]];
171    for (int u = 0; u < dim[j]; ++u) {
172      parameters[j][u] = 2.0 * RandDouble() - 1.0;
173    }
174  }
175
176  // Make a term and probe it.
177  GoodTestTerm good_term(arity, dim);
178  typedef GradientChecker<GoodTestTerm, 1, 2, 3, 4> GoodTermGradientChecker;
179  EXPECT_TRUE(GoodTermGradientChecker::Probe(
180      parameters.get(), 1e-6, &good_term, NULL));
181
182  BadTestTerm bad_term(arity, dim);
183  typedef GradientChecker<BadTestTerm, 1, 2, 3, 4> BadTermGradientChecker;
184  EXPECT_FALSE(BadTermGradientChecker::Probe(
185      parameters.get(), 1e-6, &bad_term, NULL));
186
187  for (int j = 0; j < arity; j++) {
188    delete[] parameters[j];
189  }
190}
191
192}  // namespace internal
193}  // namespace ceres
194