Change NumericDiffCostFunction to accept variadic functors.
The interface for NumericDiffCostFunction and AutoDiffCostFunction
are not comparable. They both accept variadic functors.
The change is backward compatible, as it still supports numeric
differentiation of CostFunction objects.
Some refactoring of documentation and code in auto_diff_cost_function
and its relatives was also done to make things consistent.
Change-Id: Ib5f230a1d4a85738eb187803b9c1cd7166bb3b92
diff --git a/include/ceres/autodiff_cost_function.h b/include/ceres/autodiff_cost_function.h
index bb08d64..2f86f56 100644
--- a/include/ceres/autodiff_cost_function.h
+++ b/include/ceres/autodiff_cost_function.h
@@ -28,10 +28,10 @@
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
//
-// Helpers for making CostFunctions as needed by the least squares framework,
-// with Jacobians computed via automatic differentiation. For more information
-// on automatic differentation, see the wikipedia article at
-// http://en.wikipedia.org/wiki/Automatic_differentiation
+// Create CostFunctions as needed by the least squares framework, with
+// Jacobians computed via automatic differentiation. For more
+// information on automatic differentation, see the wikipedia article
+// at http://en.wikipedia.org/wiki/Automatic_differentiation
//
// To get an auto differentiated cost function, you must define a class with a
// templated operator() (a functor) that computes the cost function in terms of
@@ -57,8 +57,8 @@
// To write an auto-differentiable cost function for the above model, first
// define the object
//
-// class MyScalarCostFunction {
-// MyScalarCostFunction(double k): k_(k) {}
+// class MyScalarCostFunctor {
+// MyScalarCostFunctor(double k): k_(k) {}
//
// template <typename T>
// bool operator()(const T* const x , const T* const y, T* e) const {
@@ -80,32 +80,32 @@
// it can be constructed as follows.
//
// CostFunction* cost_function
-// = new AutoDiffCostFunction<MyScalarCostFunction, 1, 2, 2>(
-// new MyScalarCostFunction(1.0)); ^ ^ ^
-// | | |
-// Dimension of residual ------+ | |
-// Dimension of x ----------------+ |
-// Dimension of y -------------------+
+// = new AutoDiffCostFunction<MyScalarCostFunctor, 1, 2, 2>(
+// new MyScalarCostFunctor(1.0)); ^ ^ ^
+// | | |
+// Dimension of residual -----+ | |
+// Dimension of x ---------------+ |
+// Dimension of y ------------------+
//
// In this example, there is usually an instance for each measumerent of k.
//
// In the instantiation above, the template parameters following
-// "MyScalarCostFunction", "1, 2, 2", describe the functor as computing a
+// "MyScalarCostFunctor", "1, 2, 2", describe the functor as computing a
// 1-dimensional output from two arguments, both 2-dimensional.
//
// The autodiff cost function also supports cost functions with a
// runtime-determined number of residuals. For example:
//
// CostFunction* cost_function
-// = new AutoDiffCostFunction<MyScalarCostFunction, DYNAMIC, 2, 2>(
-// new CostFunctionWithDynamicNumResiduals(1.0), ^ ^ ^
-// runtime_number_of_residuals); <----+ | | |
-// | | | |
-// | | | |
-// Actual number of residuals ------+ | | |
-// Indicate dynamic number of residuals ---------+ | |
-// Dimension of x -------------------------------------+ |
-// Dimension of y ----------------------------------------+
+// = new AutoDiffCostFunction<MyScalarCostFunctor, DYNAMIC, 2, 2>(
+// new CostFunctorWithDynamicNumResiduals(1.0), ^ ^ ^
+// runtime_number_of_residuals); <----+ | | |
+// | | | |
+// | | | |
+// Actual number of residuals ------+ | | |
+// Indicate dynamic number of residuals --------+ | |
+// Dimension of x ------------------------------------+ |
+// Dimension of y ---------------------------------------+
//
// The framework can currently accommodate cost functions of up to 6 independent
// variables, and there is no limit on the dimensionality of each of them.
@@ -119,7 +119,7 @@
// functions is to get the sizing wrong. In particular, there is a tendency to
// set the template parameters to (dimension of residual, number of parameters)
// instead of passing a dimension parameter for *every parameter*. In the
-// example above, that would be <MyScalarCostFunction, 1, 2>, which is missing
+// example above, that would be <MyScalarCostFunctor, 1, 2>, which is missing
// the last '2' argument. Please be careful when setting the size parameters.
#ifndef CERES_PUBLIC_AUTODIFF_COST_FUNCTION_H_
diff --git a/include/ceres/internal/autodiff.h b/include/ceres/internal/autodiff.h
index 36e45ac..abc16e6 100644
--- a/include/ceres/internal/autodiff.h
+++ b/include/ceres/internal/autodiff.h
@@ -146,6 +146,7 @@
#include "ceres/jet.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/fixed_array.h"
+#include "ceres/internal/variadic_evaluate.h"
namespace ceres {
namespace internal {
@@ -195,144 +196,6 @@
}
}
-// This block of quasi-repeated code calls the user-supplied functor, which may
-// take a variable number of arguments. This is accomplished by specializing the
-// struct based on the size of the trailing parameters; parameters with 0 size
-// are assumed missing.
-//
-// Supporting variadic functions is the primary source of complexity in the
-// autodiff implementation.
-
-template<typename Functor, typename T, int N0, int N1, int N2, int N3, int N4,
- int N5, int N6, int N7, int N8, int N9>
-struct VariadicEvaluate {
- static bool Call(const Functor& functor, T const *const *input, T* output) {
- return functor(input[0],
- input[1],
- input[2],
- input[3],
- input[4],
- input[5],
- input[6],
- input[7],
- input[8],
- input[9],
- output);
- }
-};
-
-template<typename Functor, typename T, int N0, int N1, int N2, int N3, int N4,
- int N5, int N6, int N7, int N8>
-struct VariadicEvaluate<Functor, T, N0, N1, N2, N3, N4, N5, N6, N7, N8, 0> {
- static bool Call(const Functor& functor, T const *const *input, T* output) {
- return functor(input[0],
- input[1],
- input[2],
- input[3],
- input[4],
- input[5],
- input[6],
- input[7],
- input[8],
- output);
- }
-};
-
-template<typename Functor, typename T, int N0, int N1, int N2, int N3, int N4,
- int N5, int N6, int N7>
-struct VariadicEvaluate<Functor, T, N0, N1, N2, N3, N4, N5, N6, N7, 0, 0> {
- static bool Call(const Functor& functor, T const *const *input, T* output) {
- return functor(input[0],
- input[1],
- input[2],
- input[3],
- input[4],
- input[5],
- input[6],
- input[7],
- output);
- }
-};
-
-template<typename Functor, typename T, int N0, int N1, int N2, int N3, int N4,
- int N5, int N6>
-struct VariadicEvaluate<Functor, T, N0, N1, N2, N3, N4, N5, N6, 0, 0, 0> {
- static bool Call(const Functor& functor, T const *const *input, T* output) {
- return functor(input[0],
- input[1],
- input[2],
- input[3],
- input[4],
- input[5],
- input[6],
- output);
- }
-};
-
-template<typename Functor, typename T, int N0, int N1, int N2, int N3, int N4,
- int N5>
-struct VariadicEvaluate<Functor, T, N0, N1, N2, N3, N4, N5, 0, 0, 0, 0> {
- static bool Call(const Functor& functor, T const *const *input, T* output) {
- return functor(input[0],
- input[1],
- input[2],
- input[3],
- input[4],
- input[5],
- output);
- }
-};
-
-template<typename Functor, typename T, int N0, int N1, int N2, int N3, int N4>
-struct VariadicEvaluate<Functor, T, N0, N1, N2, N3, N4, 0, 0, 0, 0, 0> {
- static bool Call(const Functor& functor, T const *const *input, T* output) {
- return functor(input[0],
- input[1],
- input[2],
- input[3],
- input[4],
- output);
- }
-};
-
-template<typename Functor, typename T, int N0, int N1, int N2, int N3>
-struct VariadicEvaluate<Functor, T, N0, N1, N2, N3, 0, 0, 0, 0, 0, 0> {
- static bool Call(const Functor& functor, T const *const *input, T* output) {
- return functor(input[0],
- input[1],
- input[2],
- input[3],
- output);
- }
-};
-
-template<typename Functor, typename T, int N0, int N1, int N2>
-struct VariadicEvaluate<Functor, T, N0, N1, N2, 0, 0, 0, 0, 0, 0, 0> {
- static bool Call(const Functor& functor, T const *const *input, T* output) {
- return functor(input[0],
- input[1],
- input[2],
- output);
- }
-};
-
-template<typename Functor, typename T, int N0, int N1>
-struct VariadicEvaluate<Functor, T, N0, N1, 0, 0, 0, 0, 0, 0, 0, 0> {
- static bool Call(const Functor& functor, T const *const *input, T* output) {
- return functor(input[0],
- input[1],
- output);
- }
-};
-
-template<typename Functor, typename T, int N0>
-struct VariadicEvaluate<Functor, T, N0, 0, 0, 0, 0, 0, 0, 0, 0, 0> {
- static bool Call(const Functor& functor, T const *const *input, T* output) {
- return functor(input[0],
- output);
- }
-};
-
// This is in a struct because default template parameters on a function are not
// supported in C++03 (though it is available in C++0x). N0 through N5 are the
// dimension of the input arguments to the user supplied functor.
@@ -390,6 +253,7 @@
x.get() + jet8,
x.get() + jet9,
};
+
JetT* output = x.get() + N0 + N1 + N2 + N3 + N4 + N5 + N6 + N7 + N8 + N9;
#define CERES_MAKE_1ST_ORDER_PERTURBATION(i) \
diff --git a/include/ceres/internal/numeric_diff.h b/include/ceres/internal/numeric_diff.h
new file mode 100644
index 0000000..1837df4
--- /dev/null
+++ b/include/ceres/internal/numeric_diff.h
@@ -0,0 +1,198 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2013 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+// mierle@gmail.com (Keir Mierle)
+//
+// Finite differencing routine used by NumericDiffCostFunction.
+
+#ifndef CERES_PUBLIC_INTERNAL_NUMERIC_DIFF_H_
+#define CERES_PUBLIC_INTERNAL_NUMERIC_DIFF_H_
+
+#include <cstring>
+#include <glog/logging.h>
+#include "Eigen/Dense"
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/cost_function.h"
+#include "ceres/internal/variadic_evaluate.h"
+#include "ceres/types.h"
+#include "ceres/cost_function.h"
+
+namespace ceres {
+namespace internal {
+
+// Helper templates that allow evaluation of a variadic functor or a
+// CostFunction object.
+template <typename CostFunctor,
+ int N0, int N1, int N2, int N3, int N4,
+ int N5, int N6, int N7, int N8, int N9 >
+bool EvaluateImpl(const CostFunctor* functor,
+ double const* const* parameters,
+ double* residuals,
+ const void* /* NOT USED */) {
+ return VariadicEvaluate<CostFunctor,
+ double,
+ N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>::Call(
+ *functor,
+ parameters,
+ residuals);
+}
+
+template <typename CostFunctor,
+ int N0, int N1, int N2, int N3, int N4,
+ int N5, int N6, int N7, int N8, int N9 >
+bool EvaluateImpl(const CostFunctor* functor,
+ double const* const* parameters,
+ double* residuals,
+ const CostFunction* /* NOT USED */) {
+ return functor->Evaluate(parameters, residuals, NULL);
+}
+
+// This is split from the main class because C++ doesn't allow partial template
+// specializations for member functions. The alternative is to repeat the main
+// class for differing numbers of parameters, which is also unfortunate.
+template <typename CostFunctor,
+ NumericDiffMethod kMethod,
+ int kNumResiduals,
+ int N0, int N1, int N2, int N3, int N4,
+ int N5, int N6, int N7, int N8, int N9,
+ int kParameterBlock,
+ int kParameterBlockSize>
+struct NumericDiff {
+ // Mutates parameters but must restore them before return.
+ static bool EvaluateJacobianForParameterBlock(
+ const CostFunctor* functor,
+ double const* residuals_at_eval_point,
+ const double relative_step_size,
+ double **parameters,
+ double *jacobian) {
+ using Eigen::Map;
+ using Eigen::Matrix;
+ using Eigen::RowMajor;
+ using Eigen::ColMajor;
+
+ typedef Matrix<double, kNumResiduals, 1> ResidualVector;
+ typedef Matrix<double, kParameterBlockSize, 1> ParameterVector;
+ typedef Matrix<double, kNumResiduals, kParameterBlockSize,
+ (kParameterBlockSize == 1 &&
+ kNumResiduals > 1) ? ColMajor : RowMajor> JacobianMatrix;
+
+
+ Map<JacobianMatrix> parameter_jacobian(jacobian,
+ kNumResiduals,
+ kParameterBlockSize);
+
+ // Mutate 1 element at a time and then restore.
+ Map<ParameterVector> x_plus_delta(parameters[kParameterBlock],
+ kParameterBlockSize);
+ ParameterVector x(x_plus_delta);
+ ParameterVector step_size = x.array().abs() * relative_step_size;
+
+ // To handle cases where a parameter is exactly zero, instead use
+ // the mean step_size for the other dimensions. If all the
+ // parameters are zero, there's no good answer. Take
+ // relative_step_size as a guess and hope for the best.
+ const double fallback_step_size =
+ (step_size.sum() == 0)
+ ? relative_step_size
+ : step_size.sum() / step_size.rows();
+
+ // For each parameter in the parameter block, use finite differences to
+ // compute the derivative for that parameter.
+ for (int j = 0; j < kParameterBlockSize; ++j) {
+ const double delta =
+ (step_size(j) == 0.0) ? fallback_step_size : step_size(j);
+
+ x_plus_delta(j) = x(j) + delta;
+
+ double residuals[kNumResiduals]; // NOLINT
+
+ if (!EvaluateImpl<CostFunctor, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>(
+ functor, parameters, residuals, functor)) {
+ return false;
+ }
+
+ // Compute this column of the jacobian in 3 steps:
+ // 1. Store residuals for the forward part.
+ // 2. Subtract residuals for the backward (or 0) part.
+ // 3. Divide out the run.
+ parameter_jacobian.col(j) =
+ Map<const ResidualVector>(residuals, kNumResiduals);
+
+ double one_over_delta = 1.0 / delta;
+ if (kMethod == CENTRAL) {
+ // Compute the function on the other side of x(j).
+ x_plus_delta(j) = x(j) - delta;
+
+ if (!EvaluateImpl<CostFunctor, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>(
+ functor, parameters, residuals, functor)) {
+ return false;
+ }
+
+ parameter_jacobian.col(j) -=
+ Map<ResidualVector>(residuals, kNumResiduals, 1);
+ one_over_delta /= 2;
+ } else {
+ // Forward difference only; reuse existing residuals evaluation.
+ parameter_jacobian.col(j) -=
+ Map<const ResidualVector>(residuals_at_eval_point, kNumResiduals);
+ }
+ x_plus_delta(j) = x(j); // Restore x_plus_delta.
+
+ // Divide out the run to get slope.
+ parameter_jacobian.col(j) *= one_over_delta;
+ }
+ return true;
+ }
+};
+
+template <typename CostFunctor,
+ NumericDiffMethod kMethod,
+ int kNumResiduals,
+ int N0, int N1, int N2, int N3, int N4,
+ int N5, int N6, int N7, int N8, int N9,
+ int kParameterBlock>
+struct NumericDiff<CostFunctor, kMethod, kNumResiduals,
+ N0, N1, N2, N3, N4, N5, N6, N7, N8, N9,
+ kParameterBlock, 0> {
+ // Mutates parameters but must restore them before return.
+ static bool EvaluateJacobianForParameterBlock(
+ const CostFunctor* functor,
+ double const* residuals_at_eval_point,
+ const double relative_step_size,
+ double **parameters,
+ double *jacobian) {
+ LOG(FATAL) << "Control should never reach here.";
+ return true;
+ }
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_PUBLIC_INTERNAL_NUMERIC_DIFF_H_
diff --git a/include/ceres/internal/variadic_evaluate.h b/include/ceres/internal/variadic_evaluate.h
new file mode 100644
index 0000000..4b1e4bd
--- /dev/null
+++ b/include/ceres/internal/variadic_evaluate.h
@@ -0,0 +1,182 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2013 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+// mierle@gmail.com (Keir Mierle)
+
+#ifndef CERES_PUBLIC_INTERNAL_VARIADIC_EVALUATE_H_
+#define CERES_PUBLIC_INTERNAL_VARIADIC_EVALUATE_H_
+
+#include <stddef.h>
+
+#include <glog/logging.h>
+#include "ceres/jet.h"
+#include "ceres/internal/eigen.h"
+#include "ceres/internal/fixed_array.h"
+
+namespace ceres {
+namespace internal {
+
+// This block of quasi-repeated code calls the user-supplied functor, which may
+// take a variable number of arguments. This is accomplished by specializing the
+// struct based on the size of the trailing parameters; parameters with 0 size
+// are assumed missing.
+template<typename Functor, typename T, int N0, int N1, int N2, int N3, int N4,
+ int N5, int N6, int N7, int N8, int N9>
+struct VariadicEvaluate {
+ static bool Call(const Functor& functor, T const *const *input, T* output) {
+ return functor(input[0],
+ input[1],
+ input[2],
+ input[3],
+ input[4],
+ input[5],
+ input[6],
+ input[7],
+ input[8],
+ input[9],
+ output);
+ }
+};
+
+template<typename Functor, typename T, int N0, int N1, int N2, int N3, int N4,
+ int N5, int N6, int N7, int N8>
+struct VariadicEvaluate<Functor, T, N0, N1, N2, N3, N4, N5, N6, N7, N8, 0> {
+ static bool Call(const Functor& functor, T const *const *input, T* output) {
+ return functor(input[0],
+ input[1],
+ input[2],
+ input[3],
+ input[4],
+ input[5],
+ input[6],
+ input[7],
+ input[8],
+ output);
+ }
+};
+
+template<typename Functor, typename T, int N0, int N1, int N2, int N3, int N4,
+ int N5, int N6, int N7>
+struct VariadicEvaluate<Functor, T, N0, N1, N2, N3, N4, N5, N6, N7, 0, 0> {
+ static bool Call(const Functor& functor, T const *const *input, T* output) {
+ return functor(input[0],
+ input[1],
+ input[2],
+ input[3],
+ input[4],
+ input[5],
+ input[6],
+ input[7],
+ output);
+ }
+};
+
+template<typename Functor, typename T, int N0, int N1, int N2, int N3, int N4,
+ int N5, int N6>
+struct VariadicEvaluate<Functor, T, N0, N1, N2, N3, N4, N5, N6, 0, 0, 0> {
+ static bool Call(const Functor& functor, T const *const *input, T* output) {
+ return functor(input[0],
+ input[1],
+ input[2],
+ input[3],
+ input[4],
+ input[5],
+ input[6],
+ output);
+ }
+};
+
+template<typename Functor, typename T, int N0, int N1, int N2, int N3, int N4,
+ int N5>
+struct VariadicEvaluate<Functor, T, N0, N1, N2, N3, N4, N5, 0, 0, 0, 0> {
+ static bool Call(const Functor& functor, T const *const *input, T* output) {
+ return functor(input[0],
+ input[1],
+ input[2],
+ input[3],
+ input[4],
+ input[5],
+ output);
+ }
+};
+
+template<typename Functor, typename T, int N0, int N1, int N2, int N3, int N4>
+struct VariadicEvaluate<Functor, T, N0, N1, N2, N3, N4, 0, 0, 0, 0, 0> {
+ static bool Call(const Functor& functor, T const *const *input, T* output) {
+ return functor(input[0],
+ input[1],
+ input[2],
+ input[3],
+ input[4],
+ output);
+ }
+};
+
+template<typename Functor, typename T, int N0, int N1, int N2, int N3>
+struct VariadicEvaluate<Functor, T, N0, N1, N2, N3, 0, 0, 0, 0, 0, 0> {
+ static bool Call(const Functor& functor, T const *const *input, T* output) {
+ return functor(input[0],
+ input[1],
+ input[2],
+ input[3],
+ output);
+ }
+};
+
+template<typename Functor, typename T, int N0, int N1, int N2>
+struct VariadicEvaluate<Functor, T, N0, N1, N2, 0, 0, 0, 0, 0, 0, 0> {
+ static bool Call(const Functor& functor, T const *const *input, T* output) {
+ return functor(input[0],
+ input[1],
+ input[2],
+ output);
+ }
+};
+
+template<typename Functor, typename T, int N0, int N1>
+struct VariadicEvaluate<Functor, T, N0, N1, 0, 0, 0, 0, 0, 0, 0, 0> {
+ static bool Call(const Functor& functor, T const *const *input, T* output) {
+ return functor(input[0],
+ input[1],
+ output);
+ }
+};
+
+template<typename Functor, typename T, int N0>
+struct VariadicEvaluate<Functor, T, N0, 0, 0, 0, 0, 0, 0, 0, 0, 0> {
+ static bool Call(const Functor& functor, T const *const *input, T* output) {
+ return functor(input[0],
+ output);
+ }
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_PUBLIC_INTERNAL_VARIADIC_EVALUATE_H_
diff --git a/include/ceres/numeric_diff_cost_function.h b/include/ceres/numeric_diff_cost_function.h
index 8544e44..54ba26f 100644
--- a/include/ceres/numeric_diff_cost_function.h
+++ b/include/ceres/numeric_diff_cost_function.h
@@ -27,19 +27,109 @@
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: keir@google.com (Keir Mierle)
+// sameeragarwal@google.com (Sameer Agarwal)
//
// Create CostFunctions as needed by the least squares framework with jacobians
// computed via numeric (a.k.a. finite) differentiation. For more details see
// http://en.wikipedia.org/wiki/Numerical_differentiation.
//
-// To get a numerically differentiated cost function, define a subclass of
-// CostFunction such that the Evaluate() function ignores the jacobian
-// parameter. The numeric differentiation wrapper will fill in the jacobian
-// parameter if nececssary by repeatedly calling the Evaluate() function with
-// small changes to the appropriate parameters, and computing the slope. For
-// performance, the numeric differentiation wrapper class is templated on the
-// concrete cost function, even though it could be implemented only in terms of
-// the virtual CostFunction interface.
+// To get an numerically differentiated cost function, you must define
+// a class with a operator() (a functor) that computes the residuals.
+//
+// The function must write the computed value in the last argument (the only
+// non-const one) and return true to indicate success.
+//
+// For example, consider a scalar error e = k - x'y, where both x and y are
+// two-dimensional column vector parameters, the prime sign indicates
+// transposition, and k is a constant. The form of this error, which is the
+// difference between a constant and an expression, is a common pattern in least
+// squares problems. For example, the value x'y might be the model expectation
+// for a series of measurements, where there is an instance of the cost function
+// for each measurement k.
+//
+// The actual cost added to the total problem is e^2, or (k - x'k)^2; however,
+// the squaring is implicitly done by the optimization framework.
+//
+// To write an numerically-differentiable cost function for the above model, first
+// define the object
+//
+// class MyScalarCostFunctor {
+// MyScalarCostFunctor(double k): k_(k) {}
+//
+// bool operator()(const double* const x,
+// const double* const y,
+// double* residuals) const {
+// residuals[0] = k_ - x[0] * y[0] + x[1] * y[1];
+// return true;
+// }
+//
+// private:
+// double k_;
+// };
+//
+// Note that in the declaration of operator() the input parameters x
+// and y come first, and are passed as const pointers to arrays of
+// doubles. If there were three input parameters, then the third input
+// parameter would come after y. The output is always the last
+// parameter, and is also a pointer to an array. In the example above,
+// the residual is a scalar, so only residuals[0] is set.
+//
+// Then given this class definition, the numerically differentiated
+// cost function with central differences used for computing the
+// derivative can be constructed as follows.
+//
+// CostFunction* cost_function
+// = new NumericDiffCostFunction<MyScalarCostFunctor, CENTRAL, 1, 2, 2>(
+// new MyScalarCostFunctor(1.0)); ^ ^ ^
+// | | | |
+// Finite Differencing Scheme -+ | | |
+// Dimension of residual ----------+ | |
+// Dimension of x --------------------+ |
+// Dimension of y -----------------------+
+//
+// In this example, there is usually an instance for each measumerent of k.
+//
+// In the instantiation above, the template parameters following
+// "MyScalarCostFunctor", "1, 2, 2", describe the functor as computing
+// a 1-dimensional output from two arguments, both 2-dimensional.
+//
+// The framework can currently accommodate cost functions of up to 10
+// independent variables, and there is no limit on the dimensionality
+// of each of them.
+//
+// The central difference method is considerably more accurate at the cost of
+// twice as many function evaluations than forward difference. Consider using
+// central differences begin with, and only after that works, trying forward
+// difference to improve performance.
+//
+// TODO(sameeragarwal): Add support for dynamic number of residuals.
+//
+// WARNING #1: A common beginner's error when first using
+// NumericDiffCostFunction is to get the sizing wrong. In particular,
+// there is a tendency to set the template parameters to (dimension of
+// residual, number of parameters) instead of passing a dimension
+// parameter for *every parameter*. In the example above, that would
+// be <MyScalarCostFunctor, 1, 2>, which is missing the last '2'
+// argument. Please be careful when setting the size parameters.
+//
+////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////
+//
+// ALTERNATE INTERFACE
+//
+// For a variety of reason, including compatibility with legacy code,
+// NumericDiffCostFunction can also take CostFunction objects as
+// input. The following describes how.
+//
+// To get a numerically differentiated cost function, define a
+// subclass of CostFunction such that the Evaluate() function ignores
+// the jacobian parameter. The numeric differentiation wrapper will
+// fill in the jacobian parameter if nececssary by repeatedly calling
+// the Evaluate() function with small changes to the appropriate
+// parameters, and computing the slope. For performance, the numeric
+// differentiation wrapper class is templated on the concrete cost
+// function, even though it could be implemented only in terms of the
+// virtual CostFunction interface.
//
// The numerically differentiated version of a cost function for a cost function
// can be constructed as follows:
@@ -51,233 +141,151 @@
// where MyCostFunction has 1 residual and 2 parameter blocks with sizes 4 and 8
// respectively. Look at the tests for a more detailed example.
//
-// The central difference method is considerably more accurate at the cost of
-// twice as many function evaluations than forward difference. Consider using
-// central differences begin with, and only after that works, trying forward
-// difference to improve performance.
-//
// TODO(keir): Characterize accuracy; mention pitfalls; provide alternatives.
#ifndef CERES_PUBLIC_NUMERIC_DIFF_COST_FUNCTION_H_
#define CERES_PUBLIC_NUMERIC_DIFF_COST_FUNCTION_H_
-#include <cstring>
#include <glog/logging.h>
#include "Eigen/Dense"
+#include "ceres/cost_function.h"
+#include "ceres/internal/numeric_diff.h"
#include "ceres/internal/scoped_ptr.h"
#include "ceres/sized_cost_function.h"
#include "ceres/types.h"
namespace ceres {
-enum NumericDiffMethod {
- CENTRAL,
- FORWARD
-};
-
-// This is split from the main class because C++ doesn't allow partial template
-// specializations for member functions. The alternative is to repeat the main
-// class for differing numbers of parameters, which is also unfortunate.
-template <typename CostFunctionNoJacobian,
- int num_residuals,
- int parameter_block_size,
- int parameter_block,
- NumericDiffMethod method>
-struct Differencer {
- // Mutates parameters but must restore them before return.
- static bool EvaluateJacobianForParameterBlock(
- const CostFunctionNoJacobian *function,
- double const* residuals_at_eval_point,
- double **parameters,
- double **jacobians) {
- using Eigen::Map;
- using Eigen::Matrix;
- using Eigen::RowMajor;
- using Eigen::ColMajor;
-
- typedef Matrix<double, num_residuals, 1> ResidualVector;
- typedef Matrix<double, parameter_block_size, 1> ParameterVector;
- typedef Matrix<double, num_residuals, parameter_block_size,
- (parameter_block_size == 1 &&
- num_residuals > 1) ? ColMajor : RowMajor> JacobianMatrix;
-
- Map<JacobianMatrix> parameter_jacobian(jacobians[parameter_block],
- num_residuals,
- parameter_block_size);
-
- // Mutate 1 element at a time and then restore.
- Map<ParameterVector> x_plus_delta(parameters[parameter_block],
- parameter_block_size);
- ParameterVector x(x_plus_delta);
-
- // TODO(keir): Pick a smarter number! In theory a good choice is sqrt(eps) *
- // x, which for doubles means about 1e-8 * x. However, I have found this
- // number too optimistic. This number should be exposed for users to change.
- const double kRelativeStepSize = 1e-6;
-
- ParameterVector step_size = x.array().abs() * kRelativeStepSize;
-
- // To handle cases where a parameter is exactly zero, instead use the mean
- // step_size for the other dimensions.
- double fallback_step_size = step_size.sum() / step_size.rows();
- if (fallback_step_size == 0.0) {
- // If all the parameters are zero, there's no good answer. Take
- // kRelativeStepSize as a guess and hope for the best.
- fallback_step_size = kRelativeStepSize;
- }
-
- // For each parameter in the parameter block, use finite differences to
- // compute the derivative for that parameter.
- for (int j = 0; j < parameter_block_size; ++j) {
- if (step_size(j) == 0.0) {
- // The parameter is exactly zero, so compromise and use the mean
- // step_size from the other parameters. This can break in many cases,
- // but it's hard to pick a good number without problem specific
- // knowledge.
- step_size(j) = fallback_step_size;
- }
- x_plus_delta(j) = x(j) + step_size(j);
-
- double residuals[num_residuals]; // NOLINT
- if (!function->Evaluate(parameters, residuals, NULL)) {
- // Something went wrong; bail.
- return false;
- }
-
- // Compute this column of the jacobian in 3 steps:
- // 1. Store residuals for the forward part.
- // 2. Subtract residuals for the backward (or 0) part.
- // 3. Divide out the run.
- parameter_jacobian.col(j) =
- Map<const ResidualVector>(residuals, num_residuals);
-
- double one_over_h = 1 / step_size(j);
- if (method == CENTRAL) {
- // Compute the function on the other side of x(j).
- x_plus_delta(j) = x(j) - step_size(j);
-
- if (!function->Evaluate(parameters, residuals, NULL)) {
- // Something went wrong; bail.
- return false;
- }
- parameter_jacobian.col(j) -=
- Map<ResidualVector>(residuals, num_residuals, 1);
- one_over_h /= 2;
- } else {
- // Forward difference only; reuse existing residuals evaluation.
- parameter_jacobian.col(j) -=
- Map<const ResidualVector>(residuals_at_eval_point, num_residuals);
- }
- x_plus_delta(j) = x(j); // Restore x_plus_delta.
-
- // Divide out the run to get slope.
- parameter_jacobian.col(j) *= one_over_h;
- }
- return true;
- }
-};
-
-// Prevent invalid instantiations.
-template <typename CostFunctionNoJacobian,
- int num_residuals,
- int parameter_block,
- NumericDiffMethod method>
-struct Differencer<CostFunctionNoJacobian,
- num_residuals,
- 0 /* parameter_block_size */,
- parameter_block,
- method> {
- static bool EvaluateJacobianForParameterBlock(
- const CostFunctionNoJacobian *function,
- double const* residuals_at_eval_point,
- double **parameters,
- double **jacobians) {
- LOG(FATAL) << "Shouldn't get here.";
- return true;
- }
-};
-
-template <typename CostFunctionNoJacobian,
- NumericDiffMethod method = CENTRAL, int M = 0,
- int N0 = 0, int N1 = 0, int N2 = 0, int N3 = 0, int N4 = 0, int N5 = 0>
+template <typename CostFunctor,
+ NumericDiffMethod method = CENTRAL,
+ int kNumResiduals = 0, // Number of residuals, or ceres::DYNAMIC
+ int N0 = 0, // Number of parameters in block 0.
+ int N1 = 0, // Number of parameters in block 1.
+ int N2 = 0, // Number of parameters in block 2.
+ int N3 = 0, // Number of parameters in block 3.
+ int N4 = 0, // Number of parameters in block 4.
+ int N5 = 0, // Number of parameters in block 5.
+ int N6 = 0, // Number of parameters in block 6.
+ int N7 = 0, // Number of parameters in block 7.
+ int N8 = 0, // Number of parameters in block 8.
+ int N9 = 0> // Number of parameters in block 9.
class NumericDiffCostFunction
- : public SizedCostFunction<M, N0, N1, N2, N3, N4, N5> {
+ : public SizedCostFunction<kNumResiduals, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9> {
public:
- NumericDiffCostFunction(CostFunctionNoJacobian* function,
- Ownership ownership)
- : function_(function), ownership_(ownership) {}
+ NumericDiffCostFunction(CostFunctor* functor,
+ const double relative_step_size = 1e-6)
+ :functor_(functor),
+ ownership_(TAKE_OWNERSHIP),
+ relative_step_size_(relative_step_size) {}
- virtual ~NumericDiffCostFunction() {
+ NumericDiffCostFunction(CostFunctor* functor,
+ Ownership ownership,
+ const double relative_step_size = 1e-6)
+ : functor_(functor),
+ ownership_(ownership),
+ relative_step_size_(relative_step_size) {}
+
+ ~NumericDiffCostFunction() {
if (ownership_ != TAKE_OWNERSHIP) {
- function_.release();
+ functor_.release();
}
}
virtual bool Evaluate(double const* const* parameters,
double* residuals,
double** jacobians) const {
+ using internal::FixedArray;
+ using internal::NumericDiff;
+
+ const int kNumParameters = N0 + N1 + N2 + N3 + N4 + N5 + N6 + N7 + N8 + N9;
+ const int kNumParameterBlocks =
+ (N0 > 0) + (N1 > 0) + (N2 > 0) + (N3 > 0) + (N4 > 0) +
+ (N5 > 0) + (N6 > 0) + (N7 > 0) + (N8 > 0) + (N9 > 0);
+
// Get the function value (residuals) at the the point to evaluate.
- bool success = function_->Evaluate(parameters, residuals, NULL);
- if (!success) {
- // Something went wrong; ignore the jacobian.
+ if (!internal::EvaluateImpl<CostFunctor,
+ N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>(
+ functor_.get(),
+ parameters,
+ residuals,
+ functor_.get())) {
return false;
}
+
if (!jacobians) {
- // Nothing to do; just forward.
return true;
}
// Create a copy of the parameters which will get mutated.
- const int kParametersSize = N0 + N1 + N2 + N3 + N4 + N5;
- double parameters_copy[kParametersSize];
- double *parameters_references_copy[6];
- parameters_references_copy[0] = ¶meters_copy[0];
- parameters_references_copy[1] = ¶meters_copy[0] + N0;
- parameters_references_copy[2] = ¶meters_copy[0] + N0 + N1;
- parameters_references_copy[3] = ¶meters_copy[0] + N0 + N1 + N2;
- parameters_references_copy[4] = ¶meters_copy[0] + N0 + N1 + N2 + N3;
- parameters_references_copy[5] =
- ¶meters_copy[0] + N0 + N1 + N2 + N3 + N4;
+ FixedArray<double> parameters_copy(kNumParameters);
+ FixedArray<double*> parameters_reference_copy(kNumParameterBlocks);
-#define COPY_PARAMETER_BLOCK(block) \
- if (N ## block) memcpy(parameters_references_copy[block], \
- parameters[block], \
- sizeof(double) * N ## block); // NOLINT
+ parameters_reference_copy[0] = parameters_copy.get();
+ if (N1) parameters_reference_copy[1] = parameters_reference_copy[0] + N0;
+ if (N2) parameters_reference_copy[2] = parameters_reference_copy[1] + N1;
+ if (N3) parameters_reference_copy[3] = parameters_reference_copy[2] + N2;
+ if (N4) parameters_reference_copy[4] = parameters_reference_copy[3] + N3;
+ if (N5) parameters_reference_copy[5] = parameters_reference_copy[4] + N4;
+ if (N6) parameters_reference_copy[6] = parameters_reference_copy[5] + N5;
+ if (N7) parameters_reference_copy[7] = parameters_reference_copy[6] + N6;
+ if (N7) parameters_reference_copy[8] = parameters_reference_copy[7] + N7;
+ if (N8) parameters_reference_copy[9] = parameters_reference_copy[8] + N8;
+
+#define COPY_PARAMETER_BLOCK(block) \
+ if (N ## block) memcpy(parameters_reference_copy[block], \
+ parameters[block], \
+ sizeof(double) * N ## block); // NOLINT
+
COPY_PARAMETER_BLOCK(0);
COPY_PARAMETER_BLOCK(1);
COPY_PARAMETER_BLOCK(2);
COPY_PARAMETER_BLOCK(3);
COPY_PARAMETER_BLOCK(4);
COPY_PARAMETER_BLOCK(5);
+ COPY_PARAMETER_BLOCK(6);
+ COPY_PARAMETER_BLOCK(7);
+ COPY_PARAMETER_BLOCK(8);
+ COPY_PARAMETER_BLOCK(9);
+
#undef COPY_PARAMETER_BLOCK
-#define EVALUATE_JACOBIAN_FOR_BLOCK(block) \
- if (N ## block && jacobians[block]) { \
- if (!Differencer<CostFunctionNoJacobian, /* NOLINT */ \
- M, \
- N ## block, \
- block, \
- method>::EvaluateJacobianForParameterBlock( \
- function_.get(), \
- residuals, \
- parameters_references_copy, \
- jacobians)) { \
- return false; \
- } \
+#define EVALUATE_JACOBIAN_FOR_BLOCK(block) \
+ if (N ## block && jacobians[block] != NULL) { \
+ if (!NumericDiff<CostFunctor, \
+ method, \
+ kNumResiduals, \
+ N0, N1, N2, N3, N4, N5, N6, N7, N8, N9, \
+ block, \
+ N ## block >::EvaluateJacobianForParameterBlock( \
+ functor_.get(), \
+ residuals, \
+ relative_step_size_, \
+ parameters_reference_copy.get(), \
+ jacobians[block])) { \
+ return false; \
+ } \
}
+
EVALUATE_JACOBIAN_FOR_BLOCK(0);
EVALUATE_JACOBIAN_FOR_BLOCK(1);
EVALUATE_JACOBIAN_FOR_BLOCK(2);
EVALUATE_JACOBIAN_FOR_BLOCK(3);
EVALUATE_JACOBIAN_FOR_BLOCK(4);
EVALUATE_JACOBIAN_FOR_BLOCK(5);
+ EVALUATE_JACOBIAN_FOR_BLOCK(6);
+ EVALUATE_JACOBIAN_FOR_BLOCK(7);
+ EVALUATE_JACOBIAN_FOR_BLOCK(8);
+ EVALUATE_JACOBIAN_FOR_BLOCK(9);
+
#undef EVALUATE_JACOBIAN_FOR_BLOCK
+
return true;
}
private:
- internal::scoped_ptr<CostFunctionNoJacobian> function_;
+ internal::scoped_ptr<CostFunctor> functor_;
Ownership ownership_;
+ const double relative_step_size_;
};
} // namespace ceres
diff --git a/include/ceres/types.h b/include/ceres/types.h
index edf404a..2c0f0e5 100644
--- a/include/ceres/types.h
+++ b/include/ceres/types.h
@@ -335,6 +335,11 @@
DYNAMIC = -1
};
+enum NumericDiffMethod {
+ CENTRAL,
+ FORWARD
+};
+
const char* LinearSolverTypeToString(LinearSolverType type);
bool StringToLinearSolverType(string value, LinearSolverType* type);
diff --git a/internal/ceres/numeric_diff_cost_function_test.cc b/internal/ceres/numeric_diff_cost_function_test.cc
index df12eb9..296cc22 100644
--- a/internal/ceres/numeric_diff_cost_function_test.cc
+++ b/internal/ceres/numeric_diff_cost_function_test.cc
@@ -34,7 +34,6 @@
#include <cmath>
#include <string>
#include <vector>
-#include "ceres/cost_function.h"
#include "ceres/internal/macros.h"
#include "ceres/internal/scoped_ptr.h"
#include "ceres/sized_cost_function.h"
@@ -50,48 +49,65 @@
// y1 = x1'x2 -> dy1/dx1 = x2, dy1/dx2 = x1
// y2 = (x1'x2)^2 -> dy2/dx1 = 2 * x2 * (x1'x2), dy2/dx2 = 2 * x1 * (x1'x2)
// y3 = x2'x2 -> dy3/dx1 = 0, dy3/dx2 = 2 * x2
-class TestCostFunction : public CostFunction {
- public:
- TestCostFunction() {
- set_num_residuals(3);
- mutable_parameter_block_sizes()->push_back(5); // x1.
- mutable_parameter_block_sizes()->push_back(5); // x2.
- }
- virtual bool Evaluate(double const* const* parameters,
- double* residuals,
- double** jacobians) const {
- (void) jacobians; // Ignored.
-
+struct EasyFunctor {
+ bool operator()(const double* x1, const double* x2, double* residuals) const {
residuals[0] = residuals[1] = residuals[2] = 0;
for (int i = 0; i < 5; ++i) {
- residuals[0] += parameters[0][i] * parameters[1][i];
- residuals[2] += parameters[1][i] * parameters[1][i];
+ residuals[0] += x1[i] * x2[i];
+ residuals[2] += x2[i] * x2[i];
}
residuals[1] = residuals[0] * residuals[0];
return true;
}
};
+class EasyCostFunction : public SizedCostFunction<3, 5, 5> {
+ public:
+ virtual bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) const {
+ (void) jacobians; // Ignored.
+ return EasyFunctor()(parameters[0], parameters[1], residuals);
+ }
+};
+
TEST(NumericDiffCostFunction, EasyCase) {
// Try both central and forward difference.
- internal::scoped_ptr<CostFunction> cfs[2];
+ internal::scoped_ptr<CostFunction> cfs[4];
cfs[0].reset(
- new NumericDiffCostFunction<TestCostFunction,
+ new NumericDiffCostFunction<EasyCostFunction,
CENTRAL,
3, /* number of residuals */
5, /* size of x1 */
5 /* size of x2 */>(
- new TestCostFunction, TAKE_OWNERSHIP));
+ new EasyCostFunction, TAKE_OWNERSHIP));
cfs[1].reset(
- new NumericDiffCostFunction<TestCostFunction,
+ new NumericDiffCostFunction<EasyCostFunction,
FORWARD,
3, /* number of residuals */
5, /* size of x1 */
5 /* size of x2 */>(
- new TestCostFunction, TAKE_OWNERSHIP));
+ new EasyCostFunction, TAKE_OWNERSHIP));
- for (int c = 0; c < 2; ++c) {
+ cfs[2].reset(
+ new NumericDiffCostFunction< EasyFunctor,
+ CENTRAL,
+ 3, /* number of residuals */
+ 5, /* size of x1 */
+ 5 /* size of x2 */>(
+ new EasyFunctor));
+
+ cfs[3].reset(
+ new NumericDiffCostFunction< EasyFunctor,
+ FORWARD,
+ 3, /* number of residuals */
+ 5, /* size of x1 */
+ 5 /* size of x2 */>(
+ new EasyFunctor));
+
+
+ for (int c = 0; c < 4; ++c) {
CostFunction *cost_function = cfs[c].get();
double x1[] = { 1.0, 2.0, 3.0, 4.0, 5.0 };
@@ -131,21 +147,11 @@
//
// dy1/dx1 = x2 * cos(x1'x2), dy1/dx2 = x1 * cos(x1'x2)
// dy2/dx1 = -x2 * exp(-x1'x2 / 10) / 10, dy2/dx2 = -x2 * exp(-x1'x2 / 10) / 10
-class TranscendentalTestCostFunction : public CostFunction {
- public:
- TranscendentalTestCostFunction() {
- set_num_residuals(2);
- mutable_parameter_block_sizes()->push_back(5); // x1.
- mutable_parameter_block_sizes()->push_back(5); // x2.
- }
- virtual bool Evaluate(double const* const* parameters,
- double* residuals,
- double** jacobians) const {
- (void) jacobians; // Ignored.
-
+struct TranscendentalFunctor {
+ bool operator()(const double* x1, const double* x2, double* residuals) const {
double x1x2 = 0;
for (int i = 0; i < 5; ++i) {
- x1x2 += parameters[0][i] * parameters[1][i];
+ x1x2 += x1[i] * x2[i];
}
residuals[0] = sin(x1x2);
residuals[1] = exp(-x1x2 / 10);
@@ -153,9 +159,19 @@
}
};
+class TranscendentalTestCostFunction : public SizedCostFunction<2, 5, 5> {
+ public:
+ virtual bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) const {
+ (void) jacobians; // Ignored.
+ return TranscendentalFunctor()(parameters[0], parameters[1], residuals);
+ }
+};
+
TEST(NumericDiffCostFunction, TransendentalOperationsInCostFunction) {
// Try both central and forward difference.
- internal::scoped_ptr<CostFunction> cfs[2];
+ internal::scoped_ptr<CostFunction> cfs[4];
cfs[0].reset(
new NumericDiffCostFunction<TranscendentalTestCostFunction,
CENTRAL,
@@ -172,7 +188,23 @@
5 /* size of x2 */>(
new TranscendentalTestCostFunction, TAKE_OWNERSHIP));
- for (int c = 0; c < 2; ++c) {
+ cfs[2].reset(
+ new NumericDiffCostFunction<TranscendentalFunctor,
+ CENTRAL,
+ 2, /* number of residuals */
+ 5, /* size of x1 */
+ 5 /* size of x2 */>(
+ new TranscendentalFunctor));
+
+ cfs[3].reset(
+ new NumericDiffCostFunction<TranscendentalFunctor,
+ FORWARD,
+ 2, /* number of residuals */
+ 5, /* size of x1 */
+ 5 /* size of x2 */>(
+ new TranscendentalFunctor));
+
+ for (int c = 0; c < 4; ++c) {
CostFunction *cost_function = cfs[c].get();
struct {
