Initial commit of Ceres Solver.
diff --git a/internal/ceres/runtime_numeric_diff_cost_function_test.cc b/internal/ceres/runtime_numeric_diff_cost_function_test.cc
new file mode 100644
index 0000000..6926d28
--- /dev/null
+++ b/internal/ceres/runtime_numeric_diff_cost_function_test.cc
@@ -0,0 +1,223 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2010, 2011, 2012 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: keir@google.com (Keir Mierle)
+//
+// Based on the tests in numeric_diff_cost_function.cc.
+//
+// TODO(keir): See about code duplication.
+
+#include "ceres/runtime_numeric_diff_cost_function.h"
+
+#include <algorithm>
+#include <cmath>
+#include <string>
+#include <vector>
+
+#include <glog/logging.h>
+#include "gtest/gtest.h"
+#include "ceres/stringprintf.h"
+#include "ceres/test_util.h"
+#include "ceres/cost_function.h"
+#include "ceres/internal/macros.h"
+#include "ceres/internal/scoped_ptr.h"
+
+namespace ceres {
+namespace internal {
+
+const double kRelativeEps = 1e-6;
+
+// 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.
+
+ 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[1] = residuals[0] * residuals[0];
+ return true;
+ }
+};
+
+TEST(NumericDiffCostFunction, EasyCase) {
+ // Try both central and forward difference.
+ TestCostFunction term;
+ scoped_ptr<CostFunction> cfs[2];
+ cfs[0].reset(
+ CreateRuntimeNumericDiffCostFunction(&term, CENTRAL, kRelativeEps));
+
+ cfs[1].reset(
+ CreateRuntimeNumericDiffCostFunction(&term, FORWARD, kRelativeEps));
+
+
+ for (int c = 0; c < 2; ++c) {
+ CostFunction *cost_function = cfs[c].get();
+
+ double x1[] = { 1.0, 2.0, 3.0, 4.0, 5.0 };
+ double x2[] = { 9.0, 9.0, 5.0, 5.0, 1.0 };
+ double *parameters[] = { &x1[0], &x2[0] };
+
+ double dydx1[15]; // 3 x 5, row major.
+ double dydx2[15]; // 3 x 5, row major.
+ double *jacobians[2] = { &dydx1[0], &dydx2[0] };
+
+ double residuals[3] = {-1e-100, -2e-100, -3e-100 };
+
+ ASSERT_TRUE(cost_function->Evaluate(¶meters[0],
+ &residuals[0],
+ &jacobians[0]));
+
+ EXPECT_EQ(residuals[0], 67);
+ EXPECT_EQ(residuals[1], 4489);
+ EXPECT_EQ(residuals[2], 213);
+
+ for (int i = 0; i < 5; ++i) {
+ LOG(INFO) << "c = " << c << " i = " << i;
+ const double kEps = c == 0 ? /* central */ 3e-9 : /* forward */ 2e-5;
+
+ ExpectClose(x2[i], dydx1[5 * 0 + i], kEps); // y1
+ ExpectClose(x1[i], dydx2[5 * 0 + i], kEps);
+ ExpectClose(2 * x2[i] * residuals[0], dydx1[5 * 1 + i], kEps); // y2
+ ExpectClose(2 * x1[i] * residuals[0], dydx2[5 * 1 + i], kEps);
+ ExpectClose(0.0, dydx1[5 * 2 + i], kEps); // y3
+ ExpectClose(2 * x2[i], dydx2[5 * 2 + i], kEps);
+ }
+ }
+}
+
+// y1 = sin(x1'x2)
+// y2 = exp(-x1'x2 / 10)
+//
+// 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.
+
+ double x1x2 = 0;
+ for (int i = 0; i < 5; ++i) {
+ x1x2 += parameters[0][i] * parameters[1][i];
+ }
+ residuals[0] = sin(x1x2);
+ residuals[1] = exp(-x1x2 / 10);
+ return true;
+ }
+};
+
+TEST(NumericDiffCostFunction, TransendentalOperationsInCostFunction) {
+ // Try both central and forward difference.
+ TranscendentalTestCostFunction term;
+ scoped_ptr<CostFunction> cfs[2];
+ cfs[0].reset(
+ CreateRuntimeNumericDiffCostFunction(&term, CENTRAL, kRelativeEps));
+
+ cfs[1].reset(
+ CreateRuntimeNumericDiffCostFunction(&term, FORWARD, kRelativeEps));
+
+ for (int c = 0; c < 2; ++c) {
+ CostFunction *cost_function = cfs[c].get();
+
+ struct {
+ double x1[5];
+ double x2[5];
+ } kTests[] = {
+ { { 1.0, 2.0, 3.0, 4.0, 5.0 }, // No zeros.
+ { 9.0, 9.0, 5.0, 5.0, 1.0 },
+ },
+ { { 0.0, 2.0, 3.0, 0.0, 5.0 }, // Some zeros x1.
+ { 9.0, 9.0, 5.0, 5.0, 1.0 },
+ },
+ { { 1.0, 2.0, 3.0, 1.0, 5.0 }, // Some zeros x2.
+ { 0.0, 9.0, 0.0, 5.0, 0.0 },
+ },
+ { { 0.0, 0.0, 0.0, 0.0, 0.0 }, // All zeros x1.
+ { 9.0, 9.0, 5.0, 5.0, 1.0 },
+ },
+ { { 1.0, 2.0, 3.0, 4.0, 5.0 }, // All zeros x2.
+ { 0.0, 0.0, 0.0, 0.0, 0.0 },
+ },
+ { { 0.0, 0.0, 0.0, 0.0, 0.0 }, // All zeros.
+ { 0.0, 0.0, 0.0, 0.0, 0.0 },
+ },
+ };
+ for (int k = 0; k < ARRAYSIZE(kTests); ++k) {
+ double *x1 = &(kTests[k].x1[0]);
+ double *x2 = &(kTests[k].x2[0]);
+ double *parameters[] = { x1, x2 };
+
+ double dydx1[10];
+ double dydx2[10];
+ double *jacobians[2] = { &dydx1[0], &dydx2[0] };
+
+ double residuals[2];
+
+ ASSERT_TRUE(cost_function->Evaluate(¶meters[0],
+ &residuals[0],
+ &jacobians[0]));
+ LOG(INFO) << "Ran evaluate for test k=" << k << " c=" << c;
+
+ double x1x2 = 0;
+ for (int i = 0; i < 5; ++i) {
+ x1x2 += x1[i] * x2[i];
+ }
+
+ for (int i = 0; i < 5; ++i) {
+ const double kEps = (c == 0 ? /* central */ 3e-9 : /* forward */ 2e-5);
+
+ ExpectClose( x2[i] * cos(x1x2), dydx1[5 * 0 + i], kEps); // NOLINT
+ ExpectClose( x1[i] * cos(x1x2), dydx2[5 * 0 + i], kEps); // NOLINT
+ ExpectClose(-x2[i] * exp(-x1x2 / 10.) / 10., dydx1[5 * 1 + i], kEps);
+ ExpectClose(-x1[i] * exp(-x1x2 / 10.) / 10., dydx2[5 * 1 + i], kEps);
+ }
+ }
+ }
+}
+
+} // namespace internal
+} // namespace ceres
