Squashed 'third_party/ceres/' changes from e51e9b46f..399cda773
399cda773 Update build documentation to reflect detection of Eigen via config mode
bb127272f Fix typos.
a0ec5c32a Update version history for 2.0.0RC2
3f6d27367 Unify symbol visibility configuration for all compilers
29c2912ee Unbreak the bazel build some more
bf47e1a36 Fix the Bazel build.
600e8c529 fix minor typos
bdcdcc78a update docs for changed cmake usage
3f69e5b36 Corrections from William Rucklidge
8bfdb02fb Rewrite uses of VLOG_IF and LOG_IF.
d1b35ffc1 Corrections from William Rucklidge
f34e80e91 Add dividers between licenses.
65c397dae Fix formatting
f63b1fea9 Add the MIT license text corresponding to the libmv derived files.
542613c13 minor formatting fix for trust_region_minimizer.cc
6d9e9843d Remove inclusion of ceres/eigen.h
eafeca5dc Fix a logging bug in TrustRegionMinimizer.
1fd0be916 Fix default initialisation of IterationCallback::cost
137bbe845 add info about clang-format to contributing docs
d3f66d77f fix formatting generated files (best effort)
a9c7361c8 minor formatting fix (wrongly updated in earlier commit)
7b8f675bf fix formatting for (non-generated) internal source files
921368ce3 Fix a number of typos in covariance.h
7b6b2491c fix formatting for examples
82275d8a4 some fixes for Linux and macOS install docs
9d762d74f fix formatting for public header files
c76478c48 gitignore *.pyc
4e69a475c Fix potential for mismatched release/debug TBB libraries
8e1d8e32a A number of small changes.
368a738e5 AutoDiffCostFunction: optional ownership
8cbd721c1 Add erf and erfc to jet.h, including tests in jet_test.cc
31366cff2 Benchmarks for dynamic autodiff.
29fb08aea Use CMAKE_PREFIX_PATH to pass Homebrew install location
242c703b5 Minor fixes to the documentation
79bbf9510 Add changelog for 2.0.0
41d05f13d Fix lint errors in evaluation_callback_test.cc
4b67903c1 Remove unused variables from problem_test.cc
10449fc36 Add Apache license to the LICENSE file for FixedArray
8c3ecec6d Fix some minor errors in IterationCallback docs
7d3ffcb42 Remove forced CONFIG from find_package(Eigen3)
a029fc0f9 Use latest FindTBB.cmake from VTK project
aa1abbc57 Replace use of GFLAGS_LIBRARIES with export gflags target
db2af1be8 Add Problem::EvaluateResidualBlockAssumingParametersUnchanged
ab4ed32cd Replace NULL with nullptr in the documentation.
ee280e27a Allow SubsetParameterization to accept an empty vector of constant parameters.
4b8c731d8 Fix a bug in DynamicAutoDiffCostFunction
5cb5b35a9 Fixed incorrect argument name in RotationMatrixToQuaternion()
e39d9ed1d Add a missing term and remove a superfluous word
27cab77b6 Reformulate some sentences
8ac6655ce Fix documentation formatting issues
7ef83e075 Update minimum required C++ version for Ceres to C++14
1d75e7568 Improve documentation for LocalParameterization
763398ca4 Update the section on Preconditioners
a614f788a Call EvaluationCallback before evaluating the fixed cost.
70308f7bb Simplify documentation generation.
e886d7e65 Reduce the number of minimizer iterations in evaluation_callback_test.cc
9483e6f2f Simplify DynamicCompressedRowJacobianWriter::Write
323cc55bb Update the version in package.xml to 2.0.0.
303b078b5 Fix few typos and alter a NULL to nullptr.
cca93fed6 Bypass Ceres' FindGlog.cmake in CeresConfig.cmake if possible
77fc1d0fc Use build_depend for private dependencies in Catkin package.xml
a09682f00 Fix MSVC version check to support use of clang-cl front-end
b70687fcc Add namespace qualified Ceres::ceres CMake target
99efa54bd Replace type aliases deprecated/removed in C++17/C++20 from FixedArray
adb973e4a NULL -> nullptr
27b717951 Respect FIND_QUIETLY flag in cmake config file
646959ef1 Do not export class template LineParameterization
1f128d070 Change the type of parameter index/offset to match their getter/setter
072c8f070 Initialize integer variables with integer instead of double
8c36bcc81 Use inline & -inlinehint-threshold in auto-diff benchmarks
57cf20aa5 static const -> static constexpr where we can.
40b27482a Add std::numeric_limit specialization for Jets
e751d6e4f Remove AutodiffCodegen
e9eb76f8e Remove AutodiffCodegen CMake integration
9435e08a7 More clang-tidy and wjr@ comment fixes
d93fac4b7 Remove AutodiffCodegen Tests
2281c6ed2 Fixes for comments from William Rucklidge
d797a87a4 Use Ridders' method in GradientChecker.
41675682d Fix a MSVC type deduction bug in ComputeHouseholderVector
947ec0c1f Remove AutodiffCodegen autodiff benchmarks
27183d661 Allow LocalParameterizations to have zero local size.
7ac7d79dc Remove HelloWorldCodegen example
8c8738bf8 Add photometric and relative-pose residuals to autodiff benchmarks
9f7fb66d6 Add a constant cost function to the autodiff benchmarks
ab0d373e4 Fix a comment in autodiff.h
27bb99714 Change SVD algorithm in covariance computation.
84fdac38e Add const to GetCovarianceMatrix*
6bde61d6b Add line local parameterization.
2c1c0932e Update documentation in autodiff.h
8904fa488 Inline Jet initialization in Autodiff
18a464d4e Remove an errant CR from local_parameterization.cc
5c85f2179 Use ArraySelector in Autodiff
80477ff07 Add class ArraySelector
e7a30359e Pass kNumResiduals to Autodiff
f339d71dd Refactor the automatic differentiation benchmarks.
d37b4cb15 Fix some include headers in codegen/test_utils.cc/h
550766e6d Add Autodiff Brdf Benchmark
8da9876e7 Add more autodiff benchmarks
6da364713 Fix Tukey loss function
cf4185c4e Add Codegen BA Benchmark
75dd30fae Simplify GenerateCodeForFunctor
9049688c6 Default Initialize ExpressionRef to Zero
bf1aff2f0 Fix 3+ nested Jet constructor
92d6541c7 Move Codegen files into codegen/ directory
8e962f37d Add Autodiff Codegen Tests
13c7a22ce Codegen Optimizer API
90799e29e Fix install and unnecessary string copy
032d5844c AutoDiff Code Generation - CMake Integration
d82de91b8 Add ExpressionGraph::Erase(ExpressionId)
c8e35e19f Add namespaces to generated functions and constants
75e575cae Fix use of incomplete type in defaulted Problem methods
8def19616 Remove ExpressionRef Move Constructor
f26f95410 Fix windows MSVC build.
fdf9cfd32 Add functions to find the matching ELSE, ENDIF expressions
678c05b28 Fix invert PSD matrix.
a384a7e96 Remove not used using declaration
a60136b7a Add COMMENT ExpressionType
f212c9295 Let Problem::SetParameterization be called more than once.
a3696835b use CMake function to create CeresConfigVersion
67fcff918 Make Problem movable.
19728e72d Add documentation for Problem::IsParameterBlockConstant
ba6e5fb4a Make the custom uninstall target optional
8547cbd55 Make EventLogger more efficient.
edb8322bd Update the minimum required version of Eigen to 3.3.
aa6ef417f Specify Eigen3_DIR in iOS and Android Travis CI builds
4655f2549 Use find_package() instead of find_dependency() in CeresConfig.cmake
a548766d1 Use glfags target
33dd469a5 Use Eigen3::Eigen target
47e784bb4 NULL-jacobians are handled correctly in generated autodiff code
edd54b83e Update Jet.h and rotation.h to use the new IF/ELSE macros
848c1f90c Update return type in code generator and add tests for logical functions
5010421bb Add the expression return type as a member to Expression
f4dc670ee Improve testing of the codegen system
572ec4a5a Rework Expression creation and insertion
c7337154e Disable the code generation module by default
7fa0f3db4 Explicitly state PUBLIC/PRIVATE when linking
4362a2169 Run clang-format on the public headers. Also update copyright year.
c56702aac Fix installation of codegen headers
0d03e74dc Fix the include in the autodiff codegen example
d16026440 Autodiff Codegen Part 4: Public API
d1703db45 Moved AutoDiffCodeGen macros to a separate (public) header
5ce6c063d Fix ExpressionRef copy constructor and add a move constructor
a90b5a12c Pass ExpressionRef by const reference instead of by value
ea057678c Remove MakeFunctionCall() and add test for Ternary
1084c5460 Quote all configure-expanded paths
3d756b07c Test Expressions with 'insert' instead of a macro
486d81812 Add ExpressionGraph::InsertExpression
3831a1dd3 Expression and ExpressionGraph comparison
9bb1dcb84 Remove definition of ExpressionRef::ExpressionRef(double&);
5be2e4883 Autodiff Codegen Part 3: CodeGenerator
6cd633043 Remove unused ExpressionTypes
7d0d69a4d Fix ExpressionRef
6ba8c57d2 Fix expression_test IsArithmetic
2b494cfb3 Update Travis CI to Bionic & Xcode 11.2
a3dde6877 Require Xcode >= 11.2 on macOS 10.15 (Catalina)
6fd4f072d Autodiff Codegen Part 2: Conditionals
52d6477a4 Detect and disable -fstack-check on macOS 10.15 with Xcode 11
46ca461b7 Fix `gradient_check_relative_precision` docs typo
4247d420f Autodiff Codegen Part 1: Expressions
ba62397d8 Run clang-format on jet.h
667062dcc Introduce BlockSparseMatrixData
17becf461 Remove a CHECK failure from covariance_impl.cc
d7f428e5c Add a missing cast in rotation.h
ea4d66e7e clang-tidy fixes.
be15b842a Integrate the SchurEliminatorForOneFBlock for the case <2,3,6>
087b28f1b Remove use of SetUsage as it creates compilation problems.
573046d7f Protect declarations of lapack functions under CERES_NO_LAPACK
71d638ef3 Add a specialized schur eliminator.
2ffddaccf Use override & final instead of just using virtual.
e4577dd6d Use override instead of virtual for subclasses.
3e5db5bc2 Fixing documentation typo.
82d325b73 Avoid memory allocations in Accelerate Sparse[Refactor/Solve]().
f66b51382 Fix some clang-tidy warnings.
0428e2dd0 Fix missing #include of <memory>
487c1aa51 Expose SubsetPreconditioner in the API
bf709ecac Move EvaluationCallback from Solver::Options to Problem::Options.
059bcb7f8 Drop ROS dependency on catkin
c4dbc927d Default to any other sparse libraries over Accelerate
db1f5b57a Allow some methods in Problem to use const double*.
a60c14525 Explicitly delete the copy constructor and copy assignment operator
084042c25 Lint changes from William Rucklidge
93d869020 Use selfAdjoingView<Upper> in InvertPSDMatrix.
a0cd0854a Speed up InvertPSDMatrix
7b53262b7 Allow Solver::Options::max_num_line_search_step_size_iterations = 0.
3e2cdca54 Make LineSearchMinizer work correctly with negative valued functions.
3ff12a878 Fix a clang-tidy warning in problem_test.cc
57441fe90 Fix two bugs.
1b852c57e Add Problem::EvaluateResidualBlock.
54ba6c27b Fix missing declaration warnings in Ceres code
fac46d50e Modernize ProductParameterization.
53dc6213f Add some missing string-to-enum-to-string convertors.
c0aa9a263 Add checks in rotation.h for inplace operations.
0f57fa82d Update Bazel WORKSPACE for newest Bazel
f8e5fba7b TripletSparseMatrix: guard against self-assignment
939253c20 Fix Eigen alignment issues.
bf67daf79 Add the missing <array> header to fixed_array.h
25e1cdbb6 Switch to FixedArray implementation from abseil.
d467a627b IdentityTransformation -> IdentityParameterization
eaec6a9d0 Fix more typos in CostFunctionToFunctor documentation.
99b5aa4aa Fix typos in CostFunctionToFunctor documentation.
ee7e2cb3c Set Homebrew paths via HINTS not CMAKE_PREFIX_PATH
4f8a01853 Revert "Fix custom Eigen on macos (EIGEN_INCLUDE_DIR_HINTS)"
e6c5c7226 Fix custom Eigen on macos (EIGEN_INCLUDE_DIR_HINTS)
5a56d522e Add the 3,3,3 template specialization.
df5c23116 Reorder initializer list to make -Wreorder happy
0fcfdb0b4 Fix the build breakage caused by the last commit.
9b9e9f0dc Reduce machoness of macro definition in cost_functor_to_function_test.cc
21d40daa0 Remove UTF-8 chars
9350e57a4 Enable optional use of sanitizers
0456edffb Update Travis CI Linux distro to 16.04 (Xenial)
bef0dfe35 Fix a typo in cubic_interpolation.h
056ba9bb1 Add AutoDiffFirstOrderFunction
6e527392d Update googletest/googlemock to db9b85e2.
1b2940749 Clarify documentation of BiCubicInterpolator::Evaluate for out-of-bounds values
Change-Id: Id61dd832e8fbe286deb0799aa1399d4017031dae
git-subtree-dir: third_party/ceres
git-subtree-split: 399cda773035d99eaf1f4a129a666b3c4df9d1b1
diff --git a/include/ceres/autodiff_cost_function.h b/include/ceres/autodiff_cost_function.h
index db3f6af..207f0a4 100644
--- a/include/ceres/autodiff_cost_function.h
+++ b/include/ceres/autodiff_cost_function.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -54,7 +54,7 @@
// 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 actual cost added to the total problem is e^2, or (k - x'y)^2; however,
// the squaring is implicitly done by the optimization framework.
//
// To write an auto-differentiable cost function for the above model, first
@@ -126,6 +126,7 @@
#define CERES_PUBLIC_AUTODIFF_COST_FUNCTION_H_
#include <memory>
+
#include "ceres/internal/autodiff.h"
#include "ceres/sized_cost_function.h"
#include "ceres/types.h"
@@ -152,56 +153,71 @@
int... Ns> // Number of parameters in each parameter block.
class AutoDiffCostFunction : public SizedCostFunction<kNumResiduals, Ns...> {
public:
- // Takes ownership of functor. Uses the template-provided value for the
- // number of residuals ("kNumResiduals").
- explicit AutoDiffCostFunction(CostFunctor* functor)
- : functor_(functor) {
+ // Takes ownership of functor by default. Uses the template-provided
+ // value for the number of residuals ("kNumResiduals").
+ explicit AutoDiffCostFunction(CostFunctor* functor,
+ Ownership ownership = TAKE_OWNERSHIP)
+ : functor_(functor), ownership_(ownership) {
static_assert(kNumResiduals != DYNAMIC,
"Can't run the fixed-size constructor if the number of "
"residuals is set to ceres::DYNAMIC.");
}
- // Takes ownership of functor. Ignores the template-provided
+ // Takes ownership of functor by default. Ignores the template-provided
// kNumResiduals in favor of the "num_residuals" argument provided.
//
// This allows for having autodiff cost functions which return varying
// numbers of residuals at runtime.
- AutoDiffCostFunction(CostFunctor* functor, int num_residuals)
- : functor_(functor) {
+ AutoDiffCostFunction(CostFunctor* functor,
+ int num_residuals,
+ Ownership ownership = TAKE_OWNERSHIP)
+ : functor_(functor), ownership_(ownership) {
static_assert(kNumResiduals == DYNAMIC,
"Can't run the dynamic-size constructor if the number of "
"residuals is not ceres::DYNAMIC.");
SizedCostFunction<kNumResiduals, Ns...>::set_num_residuals(num_residuals);
}
- virtual ~AutoDiffCostFunction() {}
+ explicit AutoDiffCostFunction(AutoDiffCostFunction&& other)
+ : functor_(std::move(other.functor_)), ownership_(other.ownership_) {}
+
+ virtual ~AutoDiffCostFunction() {
+ // Manually release pointer if configured to not take ownership rather than
+ // deleting only if ownership is taken.
+ // This is to stay maximally compatible to old user code which may have
+ // forgotten to implement a virtual destructor, from when the
+ // AutoDiffCostFunction always took ownership.
+ if (ownership_ == DO_NOT_TAKE_OWNERSHIP) {
+ functor_.release();
+ }
+ }
// Implementation details follow; clients of the autodiff cost function should
// not have to examine below here.
//
// To handle variadic cost functions, some template magic is needed. It's
// mostly hidden inside autodiff.h.
- virtual bool Evaluate(double const* const* parameters,
- double* residuals,
- double** jacobians) const {
+ bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) const override {
using ParameterDims =
typename SizedCostFunction<kNumResiduals, Ns...>::ParameterDims;
if (!jacobians) {
- return internal::VariadicEvaluate<ParameterDims>(*functor_,
- parameters,
- residuals);
+ return internal::VariadicEvaluate<ParameterDims>(
+ *functor_, parameters, residuals);
}
- return internal::AutoDifferentiate<ParameterDims>(
+ return internal::AutoDifferentiate<kNumResiduals, ParameterDims>(
*functor_,
parameters,
SizedCostFunction<kNumResiduals, Ns...>::num_residuals(),
residuals,
jacobians);
- }
+ };
private:
std::unique_ptr<CostFunctor> functor_;
+ Ownership ownership_;
};
} // namespace ceres
diff --git a/include/ceres/autodiff_first_order_function.h b/include/ceres/autodiff_first_order_function.h
new file mode 100644
index 0000000..b98d845
--- /dev/null
+++ b/include/ceres/autodiff_first_order_function.h
@@ -0,0 +1,151 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2019 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// 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)
+
+#ifndef CERES_PUBLIC_AUTODIFF_FIRST_ORDER_FUNCTION_H_
+#define CERES_PUBLIC_AUTODIFF_FIRST_ORDER_FUNCTION_H_
+
+#include <memory>
+
+#include "ceres/first_order_function.h"
+#include "ceres/internal/eigen.h"
+#include "ceres/internal/fixed_array.h"
+#include "ceres/jet.h"
+#include "ceres/types.h"
+
+namespace ceres {
+
+// Create FirstOrderFunctions as needed by the GradientProblem
+// framework, with gradients computed via automatic
+// differentiation. For more information on automatic differentiation,
+// see the wikipedia article at
+// http://en.wikipedia.org/wiki/Automatic_differentiation
+//
+// To get an auto differentiated function, you must define a class
+// with a templated operator() (a functor) that computes the cost
+// function in terms of the template parameter T. The autodiff
+// framework substitutes appropriate "jet" objects for T in order to
+// compute the derivative when necessary, but this is hidden, and you
+// should write the function as if T were a scalar type (e.g. a
+// double-precision floating point number).
+//
+// 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 = x'y - a, where both x and y are
+// two-dimensional column vector parameters, the prime sign indicates
+// transposition, and a is a constant.
+//
+// To write an auto-differentiable FirstOrderFunction for the above model, first
+// define the object
+//
+// class QuadraticCostFunctor {
+// public:
+// explicit QuadraticCostFunctor(double a) : a_(a) {}
+// template <typename T>
+// bool operator()(const T* const xy, T* cost) const {
+// const T* const x = xy;
+// const T* const y = xy + 2;
+// *cost = x[0] * y[0] + x[1] * y[1] - T(a_);
+// return true;
+// }
+//
+// private:
+// double a_;
+// };
+//
+// Note that in the declaration of operator() the input parameters xy come
+// first, and are passed as const pointers to arrays of T. The
+// output is the last parameter.
+//
+// Then given this class definition, the auto differentiated FirstOrderFunction
+// for it can be constructed as follows.
+//
+// FirstOrderFunction* function =
+// new AutoDiffFirstOrderFunction<QuadraticCostFunctor, 4>(
+// new QuadraticCostFunctor(1.0)));
+//
+// In the instantiation above, the template parameters following
+// "QuadraticCostFunctor", "4", describe the functor as computing a
+// 1-dimensional output from a four dimensional vector.
+//
+// WARNING: Since the functor will get instantiated with different types for
+// T, you must convert from other numeric types to T before mixing
+// computations with other variables of type T. In the example above, this is
+// seen where instead of using a_ directly, a_ is wrapped with T(a_).
+
+template <typename FirstOrderFunctor, int kNumParameters>
+class AutoDiffFirstOrderFunction : public FirstOrderFunction {
+ public:
+ // Takes ownership of functor.
+ explicit AutoDiffFirstOrderFunction(FirstOrderFunctor* functor)
+ : functor_(functor) {
+ static_assert(kNumParameters > 0, "kNumParameters must be positive");
+ }
+
+ virtual ~AutoDiffFirstOrderFunction() {}
+
+ bool Evaluate(const double* const parameters,
+ double* cost,
+ double* gradient) const override {
+ if (gradient == nullptr) {
+ return (*functor_)(parameters, cost);
+ }
+
+ typedef Jet<double, kNumParameters> JetT;
+ internal::FixedArray<JetT, (256 * 7) / sizeof(JetT)> x(kNumParameters);
+ for (int i = 0; i < kNumParameters; ++i) {
+ x[i].a = parameters[i];
+ x[i].v.setZero();
+ x[i].v[i] = 1.0;
+ }
+
+ JetT output;
+ output.a = kImpossibleValue;
+ output.v.setConstant(kImpossibleValue);
+
+ if (!(*functor_)(x.data(), &output)) {
+ return false;
+ }
+
+ *cost = output.a;
+ VectorRef(gradient, kNumParameters) = output.v;
+ return true;
+ }
+
+ int NumParameters() const override { return kNumParameters; }
+
+ private:
+ std::unique_ptr<FirstOrderFunctor> functor_;
+};
+
+} // namespace ceres
+
+#endif // CERES_PUBLIC_AUTODIFF_FIRST_ORDER_FUNCTION_H_
diff --git a/include/ceres/autodiff_local_parameterization.h b/include/ceres/autodiff_local_parameterization.h
index 649e05d..d694376 100644
--- a/include/ceres/autodiff_local_parameterization.h
+++ b/include/ceres/autodiff_local_parameterization.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -34,8 +34,9 @@
#define CERES_PUBLIC_AUTODIFF_LOCAL_PARAMETERIZATION_H_
#include <memory>
-#include "ceres/local_parameterization.h"
+
#include "ceres/internal/autodiff.h"
+#include "ceres/local_parameterization.h"
namespace ceres {
@@ -107,21 +108,20 @@
template <typename Functor, int kGlobalSize, int kLocalSize>
class AutoDiffLocalParameterization : public LocalParameterization {
public:
- AutoDiffLocalParameterization() :
- functor_(new Functor()) {}
+ AutoDiffLocalParameterization() : functor_(new Functor()) {}
// Takes ownership of functor.
- explicit AutoDiffLocalParameterization(Functor* functor) :
- functor_(functor) {}
+ explicit AutoDiffLocalParameterization(Functor* functor)
+ : functor_(functor) {}
virtual ~AutoDiffLocalParameterization() {}
- virtual bool Plus(const double* x,
- const double* delta,
- double* x_plus_delta) const {
+ bool Plus(const double* x,
+ const double* delta,
+ double* x_plus_delta) const override {
return (*functor_)(x, delta, x_plus_delta);
}
- virtual bool ComputeJacobian(const double* x, double* jacobian) const {
+ bool ComputeJacobian(const double* x, double* jacobian) const override {
double zero_delta[kLocalSize];
for (int i = 0; i < kLocalSize; ++i) {
zero_delta[i] = 0.0;
@@ -133,14 +133,15 @@
}
const double* parameter_ptrs[2] = {x, zero_delta};
- double* jacobian_ptrs[2] = { NULL, jacobian };
+ double* jacobian_ptrs[2] = {NULL, jacobian};
return internal::AutoDifferentiate<
+ kGlobalSize,
internal::StaticParameterDims<kGlobalSize, kLocalSize>>(
*functor_, parameter_ptrs, kGlobalSize, x_plus_delta, jacobian_ptrs);
}
- virtual int GlobalSize() const { return kGlobalSize; }
- virtual int LocalSize() const { return kLocalSize; }
+ int GlobalSize() const override { return kGlobalSize; }
+ int LocalSize() const override { return kLocalSize; }
private:
std::unique_ptr<Functor> functor_;
diff --git a/include/ceres/c_api.h b/include/ceres/c_api.h
index df7c9b6..91b82bf 100644
--- a/include/ceres/c_api.h
+++ b/include/ceres/c_api.h
@@ -1,5 +1,5 @@
/* Ceres Solver - A fast non-linear least squares minimizer
- * Copyright 2015 Google Inc. All rights reserved.
+ * Copyright 2019 Google Inc. All rights reserved.
* http://ceres-solver.org/
*
* Redistribution and use in source and binary forms, with or without
@@ -38,8 +38,10 @@
#ifndef CERES_PUBLIC_C_API_H_
#define CERES_PUBLIC_C_API_H_
+// clang-format off
#include "ceres/internal/port.h"
#include "ceres/internal/disable_warnings.h"
+// clang-format on
#ifdef __cplusplus
extern "C" {
@@ -143,4 +145,4 @@
#include "ceres/internal/reenable_warnings.h"
-#endif /* CERES_PUBLIC_C_API_H_ */
+#endif /* CERES_PUBLIC_C_API_H_ */
diff --git a/include/ceres/ceres.h b/include/ceres/ceres.h
index 6e077cc..d249351 100644
--- a/include/ceres/ceres.h
+++ b/include/ceres/ceres.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -46,6 +46,7 @@
#include "ceres/dynamic_cost_function.h"
#include "ceres/dynamic_cost_function_to_functor.h"
#include "ceres/dynamic_numeric_diff_cost_function.h"
+#include "ceres/evaluation_callback.h"
#include "ceres/gradient_checker.h"
#include "ceres/gradient_problem.h"
#include "ceres/gradient_problem_solver.h"
diff --git a/include/ceres/conditioned_cost_function.h b/include/ceres/conditioned_cost_function.h
index f92787e..a57ee20 100644
--- a/include/ceres/conditioned_cost_function.h
+++ b/include/ceres/conditioned_cost_function.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -34,12 +34,12 @@
#ifndef CERES_PUBLIC_CONDITIONED_COST_FUNCTION_H_
#define CERES_PUBLIC_CONDITIONED_COST_FUNCTION_H_
+#include <memory>
#include <vector>
-#include <memory>
#include "ceres/cost_function.h"
-#include "ceres/types.h"
#include "ceres/internal/disable_warnings.h"
+#include "ceres/types.h"
namespace ceres {
@@ -84,9 +84,9 @@
Ownership ownership);
virtual ~ConditionedCostFunction();
- virtual bool Evaluate(double const* const* parameters,
- double* residuals,
- double** jacobians) const;
+ bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) const override;
private:
std::unique_ptr<CostFunction> wrapped_cost_function_;
diff --git a/include/ceres/context.h b/include/ceres/context.h
index cf7c436..d08e32b 100644
--- a/include/ceres/context.h
+++ b/include/ceres/context.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2018 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
diff --git a/include/ceres/cost_function.h b/include/ceres/cost_function.h
index 39425e8..d1550c1 100644
--- a/include/ceres/cost_function.h
+++ b/include/ceres/cost_function.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2017 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -46,8 +46,9 @@
#include <cstdint>
#include <vector>
-#include "ceres/internal/port.h"
+
#include "ceres/internal/disable_warnings.h"
+#include "ceres/internal/port.h"
namespace ceres {
@@ -120,18 +121,14 @@
return parameter_block_sizes_;
}
- int num_residuals() const {
- return num_residuals_;
- }
+ int num_residuals() const { return num_residuals_; }
protected:
std::vector<int32_t>* mutable_parameter_block_sizes() {
return ¶meter_block_sizes_;
}
- void set_num_residuals(int num_residuals) {
- num_residuals_ = num_residuals;
- }
+ void set_num_residuals(int num_residuals) { num_residuals_ = num_residuals; }
private:
// Cost function signature metadata: number of inputs & their sizes,
diff --git a/include/ceres/cost_function_to_functor.h b/include/ceres/cost_function_to_functor.h
index 678c80c..9364293 100644
--- a/include/ceres/cost_function_to_functor.h
+++ b/include/ceres/cost_function_to_functor.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -38,9 +38,9 @@
// class IntrinsicProjection : public SizedCostFunction<2, 5, 3> {
// public:
// IntrinsicProjection(const double* observation);
-// virtual bool Evaluate(double const* const* parameters,
-// double* residuals,
-// double** jacobians) const;
+// bool Evaluate(double const* const* parameters,
+// double* residuals,
+// double** jacobians) const override;
// };
//
// is a cost function that implements the projection of a point in its
@@ -89,12 +89,12 @@
#include <cstdint>
#include <numeric>
#include <tuple>
+#include <utility>
#include <vector>
#include "ceres/cost_function.h"
#include "ceres/dynamic_cost_function_to_functor.h"
#include "ceres/internal/fixed_array.h"
-#include "ceres/internal/integer_sequence.h"
#include "ceres/internal/parameter_dims.h"
#include "ceres/internal/port.h"
#include "ceres/types.h"
@@ -124,8 +124,8 @@
}
}
- CHECK_EQ(accumulate(parameter_block_sizes.begin(),
- parameter_block_sizes.end(), 0),
+ CHECK_EQ(accumulate(
+ parameter_block_sizes.begin(), parameter_block_sizes.end(), 0),
ParameterDims::kNumParameters);
}
@@ -144,8 +144,7 @@
// Extract parameter block pointers from params.
using Indices =
- internal::make_integer_sequence<int,
- ParameterDims::kNumParameterBlocks>;
+ std::make_integer_sequence<int, ParameterDims::kNumParameterBlocks>;
std::array<const T*, ParameterDims::kNumParameterBlocks> parameter_blocks =
GetParameterPointers<T>(params, Indices());
@@ -158,7 +157,7 @@
template <typename T, typename Tuple, int... Indices>
static std::array<const T*, ParameterDims::kNumParameterBlocks>
GetParameterPointers(const Tuple& paramPointers,
- internal::integer_sequence<int, Indices...>) {
+ std::integer_sequence<int, Indices...>) {
return std::array<const T*, ParameterDims::kNumParameterBlocks>{
{std::get<Indices>(paramPointers)...}};
}
diff --git a/include/ceres/covariance.h b/include/ceres/covariance.h
index da9f525..2fe025d 100644
--- a/include/ceres/covariance.h
+++ b/include/ceres/covariance.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -34,6 +34,7 @@
#include <memory>
#include <utility>
#include <vector>
+
#include "ceres/internal/disable_warnings.h"
#include "ceres/internal/port.h"
#include "ceres/types.h"
@@ -50,7 +51,7 @@
// =======
// It is very easy to use this class incorrectly without understanding
// the underlying mathematics. Please read and understand the
-// documentation completely before attempting to use this class.
+// documentation completely before attempting to use it.
//
//
// This class allows the user to evaluate the covariance for a
@@ -72,7 +73,7 @@
// the maximum likelihood estimate of x given observations y is the
// solution to the non-linear least squares problem:
//
-// x* = arg min_x |f(x)|^2
+// x* = arg min_x |f(x) - y|^2
//
// And the covariance of x* is given by
//
@@ -219,11 +220,11 @@
// 1. DENSE_SVD uses Eigen's JacobiSVD to perform the
// computations. It computes the singular value decomposition
//
- // U * S * V' = J
+ // U * D * V' = J
//
// and then uses it to compute the pseudo inverse of J'J as
//
- // pseudoinverse[J'J]^ = V * pseudoinverse[S] * V'
+ // pseudoinverse[J'J] = V * pseudoinverse[D^2] * V'
//
// It is an accurate but slow method and should only be used
// for small to moderate sized problems. It can handle
@@ -234,7 +235,7 @@
//
// Q * R = J
//
- // [J'J]^-1 = [R*R']^-1
+ // [J'J]^-1 = [R'*R]^-1
//
// SPARSE_QR is not capable of computing the covariance if the
// Jacobian is rank deficient. Depending on the value of
@@ -349,10 +350,9 @@
// covariance computation. Please see the documentation for
// Covariance::Options for more on the conditions under which this
// function returns false.
- bool Compute(
- const std::vector<std::pair<const double*,
- const double*>>& covariance_blocks,
- Problem* problem);
+ bool Compute(const std::vector<std::pair<const double*, const double*>>&
+ covariance_blocks,
+ Problem* problem);
// Compute a part of the covariance matrix.
//
@@ -425,8 +425,8 @@
// a square matrix whose row and column count is equal to the sum of
// the sizes of the individual parameter blocks. The covariance
// matrix will be a row-major matrix.
- bool GetCovarianceMatrix(const std::vector<const double *> ¶meter_blocks,
- double *covariance_matrix);
+ bool GetCovarianceMatrix(const std::vector<const double*>& parameter_blocks,
+ double* covariance_matrix) const;
// Return the covariance matrix corresponding to parameter_blocks
// in the tangent space if a local parameterization is associated
@@ -445,7 +445,7 @@
// blocks. The covariance matrix will be a row-major matrix.
bool GetCovarianceMatrixInTangentSpace(
const std::vector<const double*>& parameter_blocks,
- double* covariance_matrix);
+ double* covariance_matrix) const;
private:
std::unique_ptr<internal::CovarianceImpl> impl_;
diff --git a/include/ceres/crs_matrix.h b/include/ceres/crs_matrix.h
index 23687c4..bc618fa 100644
--- a/include/ceres/crs_matrix.h
+++ b/include/ceres/crs_matrix.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -32,8 +32,9 @@
#define CERES_PUBLIC_CRS_MATRIX_H_
#include <vector>
-#include "ceres/internal/port.h"
+
#include "ceres/internal/disable_warnings.h"
+#include "ceres/internal/port.h"
namespace ceres {
diff --git a/include/ceres/cubic_interpolation.h b/include/ceres/cubic_interpolation.h
index 88a3ebf..9b9ea4a 100644
--- a/include/ceres/cubic_interpolation.h
+++ b/include/ceres/cubic_interpolation.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -31,8 +31,8 @@
#ifndef CERES_PUBLIC_CUBIC_INTERPOLATION_H_
#define CERES_PUBLIC_CUBIC_INTERPOLATION_H_
-#include "ceres/internal/port.h"
#include "Eigen/Core"
+#include "ceres/internal/port.h"
#include "glog/logging.h"
namespace ceres {
@@ -52,7 +52,7 @@
//
// "Cubic convolution interpolation for digital image processing".
// IEEE Transactions on Acoustics, Speech, and Signal Processing
-// 29 (6): 1153–1160.
+// 29 (6): 1153-1160.
//
// For more details see
//
@@ -116,15 +116,14 @@
// Example usage:
//
// const double data[] = {1.0, 2.0, 5.0, 6.0};
-// Grid1D<double, 1> grid(x, 0, 4);
+// Grid1D<double, 1> grid(data, 0, 4);
// CubicInterpolator<Grid1D<double, 1>> interpolator(grid);
// double f, dfdx;
// interpolator.Evaluator(1.5, &f, &dfdx);
-template<typename Grid>
+template <typename Grid>
class CubicInterpolator {
public:
- explicit CubicInterpolator(const Grid& grid)
- : grid_(grid) {
+ explicit CubicInterpolator(const Grid& grid) : grid_(grid) {
// The + casts the enum into an int before doing the
// comparison. It is needed to prevent
// "-Wunnamed-type-template-args" related errors.
@@ -135,7 +134,7 @@
const int n = std::floor(x);
Eigen::Matrix<double, Grid::DATA_DIMENSION, 1> p0, p1, p2, p3;
grid_.GetValue(n - 1, p0.data());
- grid_.GetValue(n, p1.data());
+ grid_.GetValue(n, p1.data());
grid_.GetValue(n + 1, p2.data());
grid_.GetValue(n + 2, p3.data());
CubicHermiteSpline<Grid::DATA_DIMENSION>(p0, p1, p2, p3, x - n, f, dfdx);
@@ -144,11 +143,10 @@
// The following two Evaluate overloads are needed for interfacing
// with automatic differentiation. The first is for when a scalar
// evaluation is done, and the second one is for when Jets are used.
- void Evaluate(const double& x, double* f) const {
- Evaluate(x, f, NULL);
- }
+ void Evaluate(const double& x, double* f) const { Evaluate(x, f, NULL); }
- template<typename JetT> void Evaluate(const JetT& x, JetT* f) const {
+ template <typename JetT>
+ void Evaluate(const JetT& x, JetT* f) const {
double fx[Grid::DATA_DIMENSION], dfdx[Grid::DATA_DIMENSION];
Evaluate(x.a, fx, dfdx);
for (int i = 0; i < Grid::DATA_DIMENSION; ++i) {
@@ -182,9 +180,7 @@
//
// f01, f11, .. fn1, f02, f12, .. , fn2
//
-template <typename T,
- int kDataDimension = 1,
- bool kInterleaved = true>
+template <typename T, int kDataDimension = 1, bool kInterleaved = true>
struct Grid1D {
public:
enum { DATA_DIMENSION = kDataDimension };
@@ -237,7 +233,7 @@
//
// "Cubic convolution interpolation for digital image processing".
// Robert G. Keys, IEEE Trans. on Acoustics, Speech, and Signal
-// Processing 29 (6): 1153–1160, 1981.
+// Processing 29 (6): 1153-1160, 1981.
//
// http://en.wikipedia.org/wiki/Cubic_Hermite_spline
// http://en.wikipedia.org/wiki/Bicubic_interpolation
@@ -252,11 +248,10 @@
// double f, dfdr, dfdc;
// interpolator.Evaluate(1.2, 2.5, &f, &dfdr, &dfdc);
-template<typename Grid>
+template <typename Grid>
class BiCubicInterpolator {
public:
- explicit BiCubicInterpolator(const Grid& grid)
- : grid_(grid) {
+ explicit BiCubicInterpolator(const Grid& grid) : grid_(grid) {
// The + casts the enum into an int before doing the
// comparison. It is needed to prevent
// "-Wunnamed-type-template-args" related errors.
@@ -264,9 +259,10 @@
}
// Evaluate the interpolated function value and/or its
- // derivative. Returns false if r or c is out of bounds.
- void Evaluate(double r, double c,
- double* f, double* dfdr, double* dfdc) const {
+ // derivative. Uses the nearest point on the grid boundary if r or
+ // c is out of bounds.
+ void Evaluate(
+ double r, double c, double* f, double* dfdr, double* dfdc) const {
// BiCubic interpolation requires 16 values around the point being
// evaluated. We will use pij, to indicate the elements of the
// 4x4 grid of values.
@@ -291,40 +287,40 @@
Eigen::Matrix<double, Grid::DATA_DIMENSION, 1> df0dc, df1dc, df2dc, df3dc;
grid_.GetValue(row - 1, col - 1, p0.data());
- grid_.GetValue(row - 1, col , p1.data());
+ grid_.GetValue(row - 1, col, p1.data());
grid_.GetValue(row - 1, col + 1, p2.data());
grid_.GetValue(row - 1, col + 2, p3.data());
- CubicHermiteSpline<Grid::DATA_DIMENSION>(p0, p1, p2, p3, c - col,
- f0.data(), df0dc.data());
+ CubicHermiteSpline<Grid::DATA_DIMENSION>(
+ p0, p1, p2, p3, c - col, f0.data(), df0dc.data());
grid_.GetValue(row, col - 1, p0.data());
- grid_.GetValue(row, col , p1.data());
+ grid_.GetValue(row, col, p1.data());
grid_.GetValue(row, col + 1, p2.data());
grid_.GetValue(row, col + 2, p3.data());
- CubicHermiteSpline<Grid::DATA_DIMENSION>(p0, p1, p2, p3, c - col,
- f1.data(), df1dc.data());
+ CubicHermiteSpline<Grid::DATA_DIMENSION>(
+ p0, p1, p2, p3, c - col, f1.data(), df1dc.data());
grid_.GetValue(row + 1, col - 1, p0.data());
- grid_.GetValue(row + 1, col , p1.data());
+ grid_.GetValue(row + 1, col, p1.data());
grid_.GetValue(row + 1, col + 1, p2.data());
grid_.GetValue(row + 1, col + 2, p3.data());
- CubicHermiteSpline<Grid::DATA_DIMENSION>(p0, p1, p2, p3, c - col,
- f2.data(), df2dc.data());
+ CubicHermiteSpline<Grid::DATA_DIMENSION>(
+ p0, p1, p2, p3, c - col, f2.data(), df2dc.data());
grid_.GetValue(row + 2, col - 1, p0.data());
- grid_.GetValue(row + 2, col , p1.data());
+ grid_.GetValue(row + 2, col, p1.data());
grid_.GetValue(row + 2, col + 1, p2.data());
grid_.GetValue(row + 2, col + 2, p3.data());
- CubicHermiteSpline<Grid::DATA_DIMENSION>(p0, p1, p2, p3, c - col,
- f3.data(), df3dc.data());
+ CubicHermiteSpline<Grid::DATA_DIMENSION>(
+ p0, p1, p2, p3, c - col, f3.data(), df3dc.data());
// Interpolate vertically the interpolated value from each row and
// compute the derivative along the columns.
CubicHermiteSpline<Grid::DATA_DIMENSION>(f0, f1, f2, f3, r - row, f, dfdr);
if (dfdc != NULL) {
// Interpolate vertically the derivative along the columns.
- CubicHermiteSpline<Grid::DATA_DIMENSION>(df0dc, df1dc, df2dc, df3dc,
- r - row, dfdc, NULL);
+ CubicHermiteSpline<Grid::DATA_DIMENSION>(
+ df0dc, df1dc, df2dc, df3dc, r - row, dfdc, NULL);
}
}
@@ -335,9 +331,8 @@
Evaluate(r, c, f, NULL, NULL);
}
- template<typename JetT> void Evaluate(const JetT& r,
- const JetT& c,
- JetT* f) const {
+ template <typename JetT>
+ void Evaluate(const JetT& r, const JetT& c, JetT* f) const {
double frc[Grid::DATA_DIMENSION];
double dfdr[Grid::DATA_DIMENSION];
double dfdc[Grid::DATA_DIMENSION];
@@ -388,12 +383,17 @@
enum { DATA_DIMENSION = kDataDimension };
Grid2D(const T* data,
- const int row_begin, const int row_end,
- const int col_begin, const int col_end)
+ const int row_begin,
+ const int row_end,
+ const int col_begin,
+ const int col_end)
: data_(data),
- row_begin_(row_begin), row_end_(row_end),
- col_begin_(col_begin), col_end_(col_end),
- num_rows_(row_end - row_begin), num_cols_(col_end - col_begin),
+ row_begin_(row_begin),
+ row_end_(row_end),
+ col_begin_(col_begin),
+ col_end_(col_end),
+ num_rows_(row_end - row_begin),
+ num_cols_(col_end - col_begin),
num_values_(num_rows_ * num_cols_) {
CHECK_GE(kDataDimension, 1);
CHECK_LT(row_begin, row_end);
@@ -406,11 +406,8 @@
const int col_idx =
std::min(std::max(col_begin_, c), col_end_ - 1) - col_begin_;
- const int n =
- (kRowMajor)
- ? num_cols_ * row_idx + col_idx
- : num_rows_ * col_idx + row_idx;
-
+ const int n = (kRowMajor) ? num_cols_ * row_idx + col_idx
+ : num_rows_ * col_idx + row_idx;
if (kInterleaved) {
for (int i = 0; i < kDataDimension; ++i) {
diff --git a/include/ceres/dynamic_autodiff_cost_function.h b/include/ceres/dynamic_autodiff_cost_function.h
index 1bfb7a5..7ccf6a8 100644
--- a/include/ceres/dynamic_autodiff_cost_function.h
+++ b/include/ceres/dynamic_autodiff_cost_function.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -33,12 +33,14 @@
#define CERES_PUBLIC_DYNAMIC_AUTODIFF_COST_FUNCTION_H_
#include <cmath>
+#include <memory>
#include <numeric>
#include <vector>
-#include <memory>
#include "ceres/dynamic_cost_function.h"
+#include "ceres/internal/fixed_array.h"
#include "ceres/jet.h"
+#include "ceres/types.h"
#include "glog/logging.h"
namespace ceres {
@@ -77,14 +79,28 @@
template <typename CostFunctor, int Stride = 4>
class DynamicAutoDiffCostFunction : public DynamicCostFunction {
public:
- explicit DynamicAutoDiffCostFunction(CostFunctor* functor)
- : functor_(functor) {}
+ // Takes ownership by default.
+ DynamicAutoDiffCostFunction(CostFunctor* functor,
+ Ownership ownership = TAKE_OWNERSHIP)
+ : functor_(functor), ownership_(ownership) {}
- virtual ~DynamicAutoDiffCostFunction() {}
+ explicit DynamicAutoDiffCostFunction(DynamicAutoDiffCostFunction&& other)
+ : functor_(std::move(other.functor_)), ownership_(other.ownership_) {}
- virtual bool Evaluate(double const* const* parameters,
- double* residuals,
- double** jacobians) const {
+ virtual ~DynamicAutoDiffCostFunction() {
+ // Manually release pointer if configured to not take ownership
+ // rather than deleting only if ownership is taken. This is to
+ // stay maximally compatible to old user code which may have
+ // forgotten to implement a virtual destructor, from when the
+ // AutoDiffCostFunction always took ownership.
+ if (ownership_ == DO_NOT_TAKE_OWNERSHIP) {
+ functor_.release();
+ }
+ }
+
+ bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) const override {
CHECK_GT(num_residuals(), 0)
<< "You must call DynamicAutoDiffCostFunction::SetNumResiduals() "
<< "before DynamicAutoDiffCostFunction::Evaluate().";
@@ -107,17 +123,19 @@
// num_residuals() derivatives. This is done with small, fixed-size jets.
const int num_parameter_blocks =
static_cast<int>(parameter_block_sizes().size());
- const int num_parameters = std::accumulate(parameter_block_sizes().begin(),
- parameter_block_sizes().end(),
- 0);
+ const int num_parameters = std::accumulate(
+ parameter_block_sizes().begin(), parameter_block_sizes().end(), 0);
// Allocate scratch space for the strided evaluation.
- std::vector<Jet<double, Stride>> input_jets(num_parameters);
- std::vector<Jet<double, Stride>> output_jets(num_residuals());
+ using JetT = Jet<double, Stride>;
+ internal::FixedArray<JetT, (256 * 7) / sizeof(JetT)> input_jets(
+ num_parameters);
+ internal::FixedArray<JetT, (256 * 7) / sizeof(JetT)> output_jets(
+ num_residuals());
// Make the parameter pack that is sent to the functor (reused).
- std::vector<Jet<double, Stride>* > jet_parameters(num_parameter_blocks,
- static_cast<Jet<double, Stride>* >(NULL));
+ internal::FixedArray<Jet<double, Stride>*> jet_parameters(
+ num_parameter_blocks, nullptr);
int num_active_parameters = 0;
// To handle constant parameters between non-constant parameter blocks, the
@@ -148,6 +166,9 @@
}
}
+ if (num_active_parameters == 0) {
+ return (*functor_)(parameters, residuals);
+ }
// When `num_active_parameters % Stride != 0` then it can be the case
// that `active_parameter_count < Stride` while parameter_cursor is less
// than the total number of parameters and with no remaining non-constant
@@ -164,8 +185,8 @@
// Evaluate all of the strides. Each stride is a chunk of the derivative to
// evaluate, typically some size proportional to the size of the SIMD
// registers of the CPU.
- int num_strides = static_cast<int>(ceil(num_active_parameters /
- static_cast<float>(Stride)));
+ int num_strides = static_cast<int>(
+ ceil(num_active_parameters / static_cast<float>(Stride)));
int current_derivative_section = 0;
int current_derivative_section_cursor = 0;
@@ -175,7 +196,7 @@
// non-constant #Stride parameters.
const int initial_derivative_section = current_derivative_section;
const int initial_derivative_section_cursor =
- current_derivative_section_cursor;
+ current_derivative_section_cursor;
int active_parameter_count = 0;
parameter_cursor = 0;
@@ -185,9 +206,9 @@
++j, parameter_cursor++) {
input_jets[parameter_cursor].v.setZero();
if (active_parameter_count < Stride &&
- parameter_cursor >= (
- start_derivative_section[current_derivative_section] +
- current_derivative_section_cursor)) {
+ parameter_cursor >=
+ (start_derivative_section[current_derivative_section] +
+ current_derivative_section_cursor)) {
if (jacobians[i] != NULL) {
input_jets[parameter_cursor].v[active_parameter_count] = 1.0;
++active_parameter_count;
@@ -214,9 +235,9 @@
for (int j = 0; j < parameter_block_sizes()[i];
++j, parameter_cursor++) {
if (active_parameter_count < Stride &&
- parameter_cursor >= (
- start_derivative_section[current_derivative_section] +
- current_derivative_section_cursor)) {
+ parameter_cursor >=
+ (start_derivative_section[current_derivative_section] +
+ current_derivative_section_cursor)) {
if (jacobians[i] != NULL) {
for (int k = 0; k < num_residuals(); ++k) {
jacobians[i][k * parameter_block_sizes()[i] + j] =
@@ -245,6 +266,7 @@
private:
std::unique_ptr<CostFunctor> functor_;
+ Ownership ownership_;
};
} // namespace ceres
diff --git a/include/ceres/dynamic_cost_function.h b/include/ceres/dynamic_cost_function.h
index 6c0aa31..6e8a076 100644
--- a/include/ceres/dynamic_cost_function.h
+++ b/include/ceres/dynamic_cost_function.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2016 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
diff --git a/include/ceres/dynamic_cost_function_to_functor.h b/include/ceres/dynamic_cost_function_to_functor.h
index 7ae5291..8d174d8 100644
--- a/include/ceres/dynamic_cost_function_to_functor.h
+++ b/include/ceres/dynamic_cost_function_to_functor.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -53,9 +53,9 @@
// class IntrinsicProjection : public CostFunction {
// public:
// IntrinsicProjection(const double* observation);
-// virtual bool Evaluate(double const* const* parameters,
-// double* residuals,
-// double** jacobians) const;
+// bool Evaluate(double const* const* parameters,
+// double* residuals,
+// double** jacobians) const override;
// };
//
// is a cost function that implements the projection of a point in its
@@ -119,8 +119,8 @@
const int num_parameter_blocks =
static_cast<int>(parameter_block_sizes.size());
const int num_residuals = cost_function_->num_residuals();
- const int num_parameters = std::accumulate(parameter_block_sizes.begin(),
- parameter_block_sizes.end(), 0);
+ const int num_parameters = std::accumulate(
+ parameter_block_sizes.begin(), parameter_block_sizes.end(), 0);
internal::FixedArray<double> parameters(num_parameters);
internal::FixedArray<double*> parameter_blocks(num_parameter_blocks);
@@ -130,8 +130,8 @@
// Build a set of arrays to get the residuals and jacobians from
// the CostFunction wrapped by this functor.
- double* parameter_ptr = parameters.get();
- double* jacobian_ptr = jacobians.get();
+ double* parameter_ptr = parameters.data();
+ double* jacobian_ptr = jacobians.data();
for (int i = 0; i < num_parameter_blocks; ++i) {
parameter_blocks[i] = parameter_ptr;
jacobian_blocks[i] = jacobian_ptr;
@@ -141,9 +141,9 @@
jacobian_ptr += num_residuals * parameter_block_sizes[i];
}
- if (!cost_function_->Evaluate(parameter_blocks.get(),
- residuals.get(),
- jacobian_blocks.get())) {
+ if (!cost_function_->Evaluate(parameter_blocks.data(),
+ residuals.data(),
+ jacobian_blocks.data())) {
return false;
}
diff --git a/include/ceres/dynamic_numeric_diff_cost_function.h b/include/ceres/dynamic_numeric_diff_cost_function.h
index 33ac5e1..ccc8f66 100644
--- a/include/ceres/dynamic_numeric_diff_cost_function.h
+++ b/include/ceres/dynamic_numeric_diff_cost_function.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -44,6 +44,7 @@
#include "ceres/internal/numeric_diff.h"
#include "ceres/internal/parameter_dims.h"
#include "ceres/numeric_diff_options.h"
+#include "ceres/types.h"
#include "glog/logging.h"
namespace ceres {
@@ -59,7 +60,9 @@
// numeric diff; the expected interface for the cost functors is:
//
// struct MyCostFunctor {
-// bool operator()(double const* const* parameters, double* residuals) const {
+// bool operator()(double const*
+// const* parameters,
+// double* residuals) const {
// // Use parameters[i] to access the i'th parameter block.
// }
// }
@@ -80,10 +83,11 @@
const CostFunctor* functor,
Ownership ownership = TAKE_OWNERSHIP,
const NumericDiffOptions& options = NumericDiffOptions())
- : functor_(functor),
- ownership_(ownership),
- options_(options) {
- }
+ : functor_(functor), ownership_(ownership), options_(options) {}
+
+ explicit DynamicNumericDiffCostFunction(
+ DynamicNumericDiffCostFunction&& other)
+ : functor_(std::move(other.functor_)), ownership_(other.ownership_) {}
virtual ~DynamicNumericDiffCostFunction() {
if (ownership_ != TAKE_OWNERSHIP) {
@@ -91,9 +95,9 @@
}
}
- virtual bool Evaluate(double const* const* parameters,
- double* residuals,
- double** jacobians) const {
+ bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) const override {
using internal::NumericDiff;
CHECK_GT(num_residuals(), 0)
<< "You must call DynamicNumericDiffCostFunction::SetNumResiduals() "
@@ -106,9 +110,7 @@
const bool status =
internal::VariadicEvaluate<internal::DynamicParameterDims>(
- *functor_.get(),
- parameters,
- residuals);
+ *functor_.get(), parameters, residuals);
if (jacobians == NULL || !status) {
return status;
}
@@ -119,8 +121,8 @@
std::vector<double*> parameters_references_copy(block_sizes.size());
parameters_references_copy[0] = ¶meters_copy[0];
for (size_t block = 1; block < block_sizes.size(); ++block) {
- parameters_references_copy[block] = parameters_references_copy[block - 1]
- + block_sizes[block - 1];
+ parameters_references_copy[block] =
+ parameters_references_copy[block - 1] + block_sizes[block - 1];
}
// Copy the parameters into the local temp space.
@@ -132,18 +134,20 @@
for (size_t block = 0; block < block_sizes.size(); ++block) {
if (jacobians[block] != NULL &&
- !NumericDiff<CostFunctor, method, ceres::DYNAMIC,
- internal::DynamicParameterDims, ceres::DYNAMIC,
+ !NumericDiff<CostFunctor,
+ method,
+ ceres::DYNAMIC,
+ internal::DynamicParameterDims,
+ ceres::DYNAMIC,
ceres::DYNAMIC>::
- EvaluateJacobianForParameterBlock(
- functor_.get(),
- residuals,
- options_,
- this->num_residuals(),
- block,
- block_sizes[block],
- ¶meters_references_copy[0],
- jacobians[block])) {
+ EvaluateJacobianForParameterBlock(functor_.get(),
+ residuals,
+ options_,
+ this->num_residuals(),
+ block,
+ block_sizes[block],
+ ¶meters_references_copy[0],
+ jacobians[block])) {
return false;
}
}
diff --git a/include/ceres/evaluation_callback.h b/include/ceres/evaluation_callback.h
index a0a8601..b9f5bbb 100644
--- a/include/ceres/evaluation_callback.h
+++ b/include/ceres/evaluation_callback.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2018 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -35,28 +35,31 @@
namespace ceres {
-// Using this callback interface, Ceres can notify you when it is about to
-// evaluate the residuals or jacobians. With the callback, you can share
-// computation between residual blocks by doing the shared computation in
-// PrepareForEvaluation() before Ceres calls CostFunction::Evaluate() on all
-// the residuals. It also enables caching results between a pure residual
-// evaluation and a residual & jacobian evaluation, via the
-// new_evaluation_point argument.
+// Using this callback interface, Ceres can notify you when it is
+// about to evaluate the residuals or jacobians. With the callback,
+// you can share computation between residual blocks by doing the
+// shared computation in PrepareForEvaluation() before Ceres calls
+// CostFunction::Evaluate(). It also enables caching results between a
+// pure residual evaluation and a residual & jacobian evaluation, via
+// the new_evaluation_point argument.
//
-// One use case for this callback is if the cost function compute is moved to
-// the GPU. In that case, the prepare call does the actual cost function
-// evaluation, and subsequent calls from Ceres to the actual cost functions
-// merely copy the results from the GPU onto the corresponding blocks for Ceres
-// to plug into the solver.
+// One use case for this callback is if the cost function compute is
+// moved to the GPU. In that case, the prepare call does the actual
+// cost function evaluation, and subsequent calls from Ceres to the
+// actual cost functions merely copy the results from the GPU onto the
+// corresponding blocks for Ceres to plug into the solver.
//
-// NOTE: Ceres provides no mechanism to share data other than the notification
-// from the callback. Users must provide access to pre-computed shared data to
-// their cost functions behind the scenes; this all happens without Ceres
-// knowing. One approach is to put a pointer to the shared data in each cost
-// function (recommended) or to use a global shared variable (discouraged;
-// bug-prone). As far as Ceres is concerned, it is evaluating cost functions
-// like any other; it just so happens that behind the scenes the cost functions
-// reuse pre-computed data to execute faster.
+// NOTE: Ceres provides no mechanism to share data other than the
+// notification from the callback. Users must provide access to
+// pre-computed shared data to their cost functions behind the scenes;
+// this all happens without Ceres knowing.
+//
+// One approach is to put a pointer to the shared data in each cost
+// function (recommended) or to use a global shared variable
+// (discouraged; bug-prone). As far as Ceres is concerned, it is
+// evaluating cost functions like any other; it just so happens that
+// behind the scenes the cost functions reuse pre-computed data to
+// execute faster.
class CERES_EXPORT EvaluationCallback {
public:
virtual ~EvaluationCallback() {}
diff --git a/include/ceres/first_order_function.h b/include/ceres/first_order_function.h
new file mode 100644
index 0000000..1420153
--- /dev/null
+++ b/include/ceres/first_order_function.h
@@ -0,0 +1,54 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2019 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// 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)
+
+#ifndef CERES_PUBLIC_FIRST_ORDER_FUNCTION_H_
+#define CERES_PUBLIC_FIRST_ORDER_FUNCTION_H_
+
+#include "ceres/internal/port.h"
+
+namespace ceres {
+
+// A FirstOrderFunction object implements the evaluation of a function
+// and its gradient.
+class CERES_EXPORT FirstOrderFunction {
+ public:
+ virtual ~FirstOrderFunction() {}
+
+ // cost is never null. gradient may be null. The return value
+ // indicates whether the evaluation was successful or not.
+ virtual bool Evaluate(const double* const parameters,
+ double* cost,
+ double* gradient) const = 0;
+ virtual int NumParameters() const = 0;
+};
+
+} // namespace ceres
+
+#endif // CERES_PUBLIC_FIRST_ORDER_FUNCTION_H_
diff --git a/include/ceres/gradient_checker.h b/include/ceres/gradient_checker.h
index e23df76..b79cf86 100644
--- a/include/ceres/gradient_checker.h
+++ b/include/ceres/gradient_checker.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
diff --git a/include/ceres/gradient_problem.h b/include/ceres/gradient_problem.h
index 6adcfd0..49d605e 100644
--- a/include/ceres/gradient_problem.h
+++ b/include/ceres/gradient_problem.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -32,6 +32,8 @@
#define CERES_PUBLIC_GRADIENT_PROBLEM_H_
#include <memory>
+
+#include "ceres/first_order_function.h"
#include "ceres/internal/port.h"
#include "ceres/local_parameterization.h"
@@ -109,18 +111,6 @@
std::unique_ptr<double[]> scratch_;
};
-// A FirstOrderFunction object implements the evaluation of a function
-// and its gradient.
-class CERES_EXPORT FirstOrderFunction {
- public:
- virtual ~FirstOrderFunction() {}
- // cost is never NULL. gradient may be null.
- virtual bool Evaluate(const double* const parameters,
- double* cost,
- double* gradient) const = 0;
- virtual int NumParameters() const = 0;
-};
-
} // namespace ceres
#endif // CERES_PUBLIC_GRADIENT_PROBLEM_H_
diff --git a/include/ceres/gradient_problem_solver.h b/include/ceres/gradient_problem_solver.h
index 1831d8d..9fab62e 100644
--- a/include/ceres/gradient_problem_solver.h
+++ b/include/ceres/gradient_problem_solver.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -34,10 +34,11 @@
#include <cmath>
#include <string>
#include <vector>
+
+#include "ceres/internal/disable_warnings.h"
#include "ceres/internal/port.h"
#include "ceres/iteration_callback.h"
#include "ceres/types.h"
-#include "ceres/internal/disable_warnings.h"
namespace ceres {
@@ -87,7 +88,7 @@
// method, please see:
//
// Nocedal, J. (1980). "Updating Quasi-Newton Matrices with
- // Limited Storage". Mathematics of Computation 35 (151): 773–782.
+ // Limited Storage". Mathematics of Computation 35 (151): 773-782.
int max_lbfgs_rank = 20;
// As part of the (L)BFGS update step (BFGS) / right-multiply step (L-BFGS),
@@ -319,7 +320,8 @@
// If the line search direction is NONLINEAR_CONJUGATE_GRADIENT,
// then this indicates the particular variant of non-linear
// conjugate gradient used.
- NonlinearConjugateGradientType nonlinear_conjugate_gradient_type = FLETCHER_REEVES;
+ NonlinearConjugateGradientType nonlinear_conjugate_gradient_type =
+ FLETCHER_REEVES;
// If the type of the line search direction is LBFGS, then this
// indicates the rank of the Hessian approximation.
diff --git a/include/ceres/internal/array_selector.h b/include/ceres/internal/array_selector.h
new file mode 100644
index 0000000..841797f
--- /dev/null
+++ b/include/ceres/internal/array_selector.h
@@ -0,0 +1,95 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2020 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// 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: darius.rueckert@fau.de (Darius Rueckert)
+//
+
+#ifndef CERES_PUBLIC_INTERNAL_ARRAY_SELECTOR_H_
+#define CERES_PUBLIC_INTERNAL_ARRAY_SELECTOR_H_
+
+#include <array>
+#include <vector>
+
+#include "ceres/internal/fixed_array.h"
+#include "ceres/types.h"
+
+namespace ceres {
+namespace internal {
+
+// StaticFixedArray selects the best array implementation based on template
+// arguments. If the size is not known at compile-time, pass
+// ceres::DYNAMIC as a size-template argument.
+//
+// Three different containers are selected in different scenarios:
+//
+// num_elements == DYNAMIC:
+// -> ceres::internal::FixedArray<T, max_stack_size>(size)
+
+// num_elements != DYNAMIC && num_elements <= max_stack_size
+// -> std::array<T,num_elements>
+
+// num_elements != DYNAMIC && num_elements > max_stack_size
+// -> std::vector<T>(num_elements)
+//
+template <typename T,
+ int num_elements,
+ int max_num_elements_on_stack,
+ bool dynamic = (num_elements == DYNAMIC),
+ bool fits_on_stack = (num_elements <= max_num_elements_on_stack)>
+struct ArraySelector {};
+
+template <typename T,
+ int num_elements,
+ int max_num_elements_on_stack,
+ bool fits_on_stack>
+struct ArraySelector<T,
+ num_elements,
+ max_num_elements_on_stack,
+ true,
+ fits_on_stack>
+ : ceres::internal::FixedArray<T, max_num_elements_on_stack> {
+ ArraySelector(int s)
+ : ceres::internal::FixedArray<T, max_num_elements_on_stack>(s) {}
+};
+
+template <typename T, int num_elements, int max_num_elements_on_stack>
+struct ArraySelector<T, num_elements, max_num_elements_on_stack, false, true>
+ : std::array<T, num_elements> {
+ ArraySelector(int s) { CHECK_EQ(s, num_elements); }
+};
+
+template <typename T, int num_elements, int max_num_elements_on_stack>
+struct ArraySelector<T, num_elements, max_num_elements_on_stack, false, false>
+ : std::vector<T> {
+ ArraySelector(int s) : std::vector<T>(s) { CHECK_EQ(s, num_elements); }
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_PUBLIC_INTERNAL_ARRAY_SELECTOR_H_
diff --git a/include/ceres/internal/autodiff.h b/include/ceres/internal/autodiff.h
index ff47fbf..9d7de75 100644
--- a/include/ceres/internal/autodiff.h
+++ b/include/ceres/internal/autodiff.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -143,7 +143,9 @@
#include <stddef.h>
#include <array>
+#include <utility>
+#include "ceres/internal/array_selector.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/fixed_array.h"
#include "ceres/internal/parameter_dims.h"
@@ -152,6 +154,17 @@
#include "ceres/types.h"
#include "glog/logging.h"
+// If the number of parameters exceeds this values, the corresponding jets are
+// placed on the heap. This will reduce performance by a factor of 2-5 on
+// current compilers.
+#ifndef CERES_AUTODIFF_MAX_PARAMETERS_ON_STACK
+#define CERES_AUTODIFF_MAX_PARAMETERS_ON_STACK 50
+#endif
+
+#ifndef CERES_AUTODIFF_MAX_RESIDUALS_ON_STACK
+#define CERES_AUTODIFF_MAX_RESIDUALS_ON_STACK 20
+#endif
+
namespace ceres {
namespace internal {
@@ -169,16 +182,24 @@
//
// is what would get put in dst if N was 3, offset was 3, and the jet type JetT
// was 8-dimensional.
-template <int Offset, int N, typename T, typename JetT>
-inline void Make1stOrderPerturbation(const T* src, JetT* dst) {
- DCHECK(src);
- DCHECK(dst);
- for (int j = 0; j < N; ++j) {
- dst[j].a = src[j];
- dst[j].v.setZero();
- dst[j].v[Offset + j] = T(1.0);
+template <int j, int N, int Offset, typename T, typename JetT>
+struct Make1stOrderPerturbation {
+ public:
+ inline static void Apply(const T* src, JetT* dst) {
+ if (j == 0) {
+ DCHECK(src);
+ DCHECK(dst);
+ }
+ dst[j] = JetT(src[j], j + Offset);
+ Make1stOrderPerturbation<j + 1, N, Offset, T, JetT>::Apply(src, dst);
}
-}
+};
+
+template <int N, int Offset, typename T, typename JetT>
+struct Make1stOrderPerturbation<N, N, Offset, T, JetT> {
+ public:
+ static void Apply(const T* src, JetT* dst) {}
+};
// Calls Make1stOrderPerturbation for every parameter block.
//
@@ -186,26 +207,31 @@
// If one having three parameter blocks with dimensions (3, 2, 4), the call
// Make1stOrderPerturbations<integer_sequence<3, 2, 4>::Apply(params, x);
// will result in the following calls to Make1stOrderPerturbation:
-// Make1stOrderPerturbation<0, 3>(params[0], x + 0);
-// Make1stOrderPerturbation<3, 2>(params[1], x + 3);
-// Make1stOrderPerturbation<5, 4>(params[2], x + 5);
+// Make1stOrderPerturbation<0, 3, 0>::Apply(params[0], x + 0);
+// Make1stOrderPerturbation<0, 2, 3>::Apply(params[1], x + 3);
+// Make1stOrderPerturbation<0, 4, 5>::Apply(params[2], x + 5);
template <typename Seq, int ParameterIdx = 0, int Offset = 0>
struct Make1stOrderPerturbations;
template <int N, int... Ns, int ParameterIdx, int Offset>
-struct Make1stOrderPerturbations<integer_sequence<int, N, Ns...>, ParameterIdx,
+struct Make1stOrderPerturbations<std::integer_sequence<int, N, Ns...>,
+ ParameterIdx,
Offset> {
template <typename T, typename JetT>
- static void Apply(T const* const* parameters, JetT* x) {
- Make1stOrderPerturbation<Offset, N>(parameters[ParameterIdx], x + Offset);
- Make1stOrderPerturbations<integer_sequence<int, Ns...>, ParameterIdx + 1,
+ inline static void Apply(T const* const* parameters, JetT* x) {
+ Make1stOrderPerturbation<0, N, Offset, T, JetT>::Apply(
+ parameters[ParameterIdx], x + Offset);
+ Make1stOrderPerturbations<std::integer_sequence<int, Ns...>,
+ ParameterIdx + 1,
Offset + N>::Apply(parameters, x);
}
};
// End of 'recursion'. Nothing more to do.
template <int ParameterIdx, int Total>
-struct Make1stOrderPerturbations<integer_sequence<int>, ParameterIdx, Total> {
+struct Make1stOrderPerturbations<std::integer_sequence<int>,
+ ParameterIdx,
+ Total> {
template <typename T, typename JetT>
static void Apply(T const* const* /* NOT USED */, JetT* /* NOT USED */) {}
};
@@ -253,61 +279,83 @@
struct Take1stOrderParts;
template <int N, int... Ns, int ParameterIdx, int Offset>
-struct Take1stOrderParts<integer_sequence<int, N, Ns...>, ParameterIdx,
+struct Take1stOrderParts<std::integer_sequence<int, N, Ns...>,
+ ParameterIdx,
Offset> {
template <typename JetT, typename T>
- static void Apply(int num_outputs, JetT* output, T** jacobians) {
+ inline static void Apply(int num_outputs, JetT* output, T** jacobians) {
if (jacobians[ParameterIdx]) {
Take1stOrderPart<Offset, N>(num_outputs, output, jacobians[ParameterIdx]);
}
- Take1stOrderParts<integer_sequence<int, Ns...>, ParameterIdx + 1,
+ Take1stOrderParts<std::integer_sequence<int, Ns...>,
+ ParameterIdx + 1,
Offset + N>::Apply(num_outputs, output, jacobians);
}
};
// End of 'recursion'. Nothing more to do.
template <int ParameterIdx, int Offset>
-struct Take1stOrderParts<integer_sequence<int>, ParameterIdx, Offset> {
+struct Take1stOrderParts<std::integer_sequence<int>, ParameterIdx, Offset> {
template <typename T, typename JetT>
- static void Apply(int /* NOT USED*/, JetT* /* NOT USED*/,
+ static void Apply(int /* NOT USED*/,
+ JetT* /* NOT USED*/,
T** /* NOT USED */) {}
};
-template <typename ParameterDims, typename Functor, typename T>
+template <int kNumResiduals,
+ typename ParameterDims,
+ typename Functor,
+ typename T>
inline bool AutoDifferentiate(const Functor& functor,
- T const *const *parameters,
- int num_outputs,
+ T const* const* parameters,
+ int dynamic_num_outputs,
T* function_value,
T** jacobians) {
- DCHECK_GT(num_outputs, 0);
-
typedef Jet<T, ParameterDims::kNumParameters> JetT;
- FixedArray<JetT, (256 * 7) / sizeof(JetT)> x(ParameterDims::kNumParameters +
- num_outputs);
-
using Parameters = typename ParameterDims::Parameters;
- // These are the positions of the respective jets in the fixed array x.
+ if (kNumResiduals != DYNAMIC) {
+ DCHECK_EQ(kNumResiduals, dynamic_num_outputs);
+ }
+
+ ArraySelector<JetT,
+ ParameterDims::kNumParameters,
+ CERES_AUTODIFF_MAX_PARAMETERS_ON_STACK>
+ parameters_as_jets(ParameterDims::kNumParameters);
+
+ // Pointers to the beginning of each parameter block
std::array<JetT*, ParameterDims::kNumParameterBlocks> unpacked_parameters =
- ParameterDims::GetUnpackedParameters(x.get());
- JetT* output = x.get() + ParameterDims::kNumParameters;
+ ParameterDims::GetUnpackedParameters(parameters_as_jets.data());
+
+ // If the number of residuals is fixed, we use the template argument as the
+ // number of outputs. Otherwise we use the num_outputs parameter. Note: The
+ // ?-operator here is compile-time evaluated, therefore num_outputs is also
+ // a compile-time constant for functors with fixed residuals.
+ const int num_outputs =
+ kNumResiduals == DYNAMIC ? dynamic_num_outputs : kNumResiduals;
+ DCHECK_GT(num_outputs, 0);
+
+ ArraySelector<JetT, kNumResiduals, CERES_AUTODIFF_MAX_RESIDUALS_ON_STACK>
+ residuals_as_jets(num_outputs);
// Invalidate the output Jets, so that we can detect if the user
// did not assign values to all of them.
for (int i = 0; i < num_outputs; ++i) {
- output[i].a = kImpossibleValue;
- output[i].v.setConstant(kImpossibleValue);
+ residuals_as_jets[i].a = kImpossibleValue;
+ residuals_as_jets[i].v.setConstant(kImpossibleValue);
}
- Make1stOrderPerturbations<Parameters>::Apply(parameters, x.get());
+ Make1stOrderPerturbations<Parameters>::Apply(parameters,
+ parameters_as_jets.data());
- if (!VariadicEvaluate<ParameterDims>(functor, unpacked_parameters.data(),
- output)) {
+ if (!VariadicEvaluate<ParameterDims>(
+ functor, unpacked_parameters.data(), residuals_as_jets.data())) {
return false;
}
- Take0thOrderPart(num_outputs, output, function_value);
- Take1stOrderParts<Parameters>::Apply(num_outputs, output, jacobians);
+ Take0thOrderPart(num_outputs, residuals_as_jets.data(), function_value);
+ Take1stOrderParts<Parameters>::Apply(
+ num_outputs, residuals_as_jets.data(), jacobians);
return true;
}
diff --git a/include/ceres/internal/disable_warnings.h b/include/ceres/internal/disable_warnings.h
index fd848fe..d7766a0 100644
--- a/include/ceres/internal/disable_warnings.h
+++ b/include/ceres/internal/disable_warnings.h
@@ -34,11 +34,11 @@
#define CERES_WARNINGS_DISABLED
#ifdef _MSC_VER
-#pragma warning( push )
+#pragma warning(push)
// Disable the warning C4251 which is triggered by stl classes in
// Ceres' public interface. To quote MSDN: "C4251 can be ignored "
// "if you are deriving from a type in the Standard C++ Library"
-#pragma warning( disable : 4251 )
+#pragma warning(disable : 4251)
#endif
#endif // CERES_WARNINGS_DISABLED
diff --git a/include/ceres/internal/eigen.h b/include/ceres/internal/eigen.h
index 59545df..b6d0b7f 100644
--- a/include/ceres/internal/eigen.h
+++ b/include/ceres/internal/eigen.h
@@ -36,31 +36,26 @@
namespace ceres {
typedef Eigen::Matrix<double, Eigen::Dynamic, 1> Vector;
-typedef Eigen::Matrix<double,
- Eigen::Dynamic,
- Eigen::Dynamic,
- Eigen::RowMajor> Matrix;
+typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
+ Matrix;
typedef Eigen::Map<Vector> VectorRef;
typedef Eigen::Map<Matrix> MatrixRef;
typedef Eigen::Map<const Vector> ConstVectorRef;
typedef Eigen::Map<const Matrix> ConstMatrixRef;
// Column major matrices for DenseSparseMatrix/DenseQRSolver
-typedef Eigen::Matrix<double,
- Eigen::Dynamic,
- Eigen::Dynamic,
- Eigen::ColMajor> ColMajorMatrix;
+typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
+ ColMajorMatrix;
-typedef Eigen::Map<ColMajorMatrix, 0,
- Eigen::Stride<Eigen::Dynamic, 1>> ColMajorMatrixRef;
+typedef Eigen::Map<ColMajorMatrix, 0, Eigen::Stride<Eigen::Dynamic, 1>>
+ ColMajorMatrixRef;
-typedef Eigen::Map<const ColMajorMatrix,
- 0,
- Eigen::Stride<Eigen::Dynamic, 1>> ConstColMajorMatrixRef;
+typedef Eigen::Map<const ColMajorMatrix, 0, Eigen::Stride<Eigen::Dynamic, 1>>
+ ConstColMajorMatrixRef;
// C++ does not support templated typdefs, thus the need for this
// struct so that we can support statically sized Matrix and Maps.
- template <int num_rows = Eigen::Dynamic, int num_cols = Eigen::Dynamic>
+template <int num_rows = Eigen::Dynamic, int num_cols = Eigen::Dynamic>
struct EigenTypes {
typedef Eigen::Matrix<double,
num_rows,
diff --git a/include/ceres/internal/fixed_array.h b/include/ceres/internal/fixed_array.h
index dcb2ace..dcbddcd 100644
--- a/include/ceres/internal/fixed_array.h
+++ b/include/ceres/internal/fixed_array.h
@@ -1,188 +1,465 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
-// http://ceres-solver.org/
+// Copyright 2018 The Abseil Authors.
//
-// Redistribution and use in source and binary forms, with or without
-// modification, are permitted provided that the following conditions are met:
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
//
-// * 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.
+// https://www.apache.org/licenses/LICENSE-2.0
//
-// 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.
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
//
-// Author: rennie@google.com (Jeffrey Rennie)
-// Author: sanjay@google.com (Sanjay Ghemawat) -- renamed to FixedArray
+// -----------------------------------------------------------------------------
+// File: fixed_array.h
+// -----------------------------------------------------------------------------
+//
+// A `FixedArray<T>` represents a non-resizable array of `T` where the length of
+// the array can be determined at run-time. It is a good replacement for
+// non-standard and deprecated uses of `alloca()` and variable length arrays
+// within the GCC extension. (See
+// https://gcc.gnu.org/onlinedocs/gcc/Variable-Length.html).
+//
+// `FixedArray` allocates small arrays inline, keeping performance fast by
+// avoiding heap operations. It also helps reduce the chances of
+// accidentally overflowing your stack if large input is passed to
+// your function.
#ifndef CERES_PUBLIC_INTERNAL_FIXED_ARRAY_H_
#define CERES_PUBLIC_INTERNAL_FIXED_ARRAY_H_
+#include <Eigen/Core> // For Eigen::aligned_allocator
+#include <algorithm>
+#include <array>
#include <cstddef>
-#include "Eigen/Core"
-#include "ceres/internal/manual_constructor.h"
+#include <memory>
+#include <tuple>
+#include <type_traits>
+
+#include "ceres/internal/memory.h"
#include "glog/logging.h"
namespace ceres {
namespace internal {
-// A FixedArray<T> represents a non-resizable array of T where the
-// length of the array does not need to be a compile time constant.
-//
-// FixedArray allocates small arrays inline, and large arrays on
-// the heap. It is a good replacement for non-standard and deprecated
-// uses of alloca() and variable length arrays (a GCC extension).
-//
-// FixedArray keeps performance fast for small arrays, because it
-// avoids heap operations. It also helps reduce the chances of
-// accidentally overflowing your stack if large input is passed to
-// your function.
-//
-// Also, FixedArray is useful for writing portable code. Not all
-// compilers support arrays of dynamic size.
+constexpr static auto kFixedArrayUseDefault = static_cast<size_t>(-1);
-// Most users should not specify an inline_elements argument and let
-// FixedArray<> automatically determine the number of elements
-// to store inline based on sizeof(T).
+// The default fixed array allocator.
//
-// If inline_elements is specified, the FixedArray<> implementation
-// will store arrays of length <= inline_elements inline.
-//
-// Finally note that unlike vector<T> FixedArray<T> will not zero-initialize
-// simple types like int, double, bool, etc.
-//
-// Non-POD types will be default-initialized just like regular vectors or
-// arrays.
+// As one can not easily detect if a struct contains or inherits from a fixed
+// size Eigen type, to be safe the Eigen::aligned_allocator is used by default.
+// But trivial types can never contain Eigen types, so std::allocator is used to
+// safe some heap memory.
+template <typename T>
+using FixedArrayDefaultAllocator =
+ typename std::conditional<std::is_trivial<T>::value,
+ std::allocator<T>,
+ Eigen::aligned_allocator<T>>::type;
-#if defined(_WIN64)
- typedef __int64 ssize_t;
-#elif defined(_WIN32)
- typedef __int32 ssize_t;
-#endif
-
-template <typename T, ssize_t inline_elements = -1>
+// -----------------------------------------------------------------------------
+// FixedArray
+// -----------------------------------------------------------------------------
+//
+// A `FixedArray` provides a run-time fixed-size array, allocating a small array
+// inline for efficiency.
+//
+// Most users should not specify an `inline_elements` argument and let
+// `FixedArray` automatically determine the number of elements
+// to store inline based on `sizeof(T)`. If `inline_elements` is specified, the
+// `FixedArray` implementation will use inline storage for arrays with a
+// length <= `inline_elements`.
+//
+// Note that a `FixedArray` constructed with a `size_type` argument will
+// default-initialize its values by leaving trivially constructible types
+// uninitialized (e.g. int, int[4], double), and others default-constructed.
+// This matches the behavior of c-style arrays and `std::array`, but not
+// `std::vector`.
+//
+// Note that `FixedArray` does not provide a public allocator; if it requires a
+// heap allocation, it will do so with global `::operator new[]()` and
+// `::operator delete[]()`, even if T provides class-scope overrides for these
+// operators.
+template <typename T,
+ size_t N = kFixedArrayUseDefault,
+ typename A = FixedArrayDefaultAllocator<T>>
class FixedArray {
+ static_assert(!std::is_array<T>::value || std::extent<T>::value > 0,
+ "Arrays with unknown bounds cannot be used with FixedArray.");
+
+ static constexpr size_t kInlineBytesDefault = 256;
+
+ using AllocatorTraits = std::allocator_traits<A>;
+ // std::iterator_traits isn't guaranteed to be SFINAE-friendly until C++17,
+ // but this seems to be mostly pedantic.
+ template <typename Iterator>
+ using EnableIfForwardIterator = typename std::enable_if<std::is_convertible<
+ typename std::iterator_traits<Iterator>::iterator_category,
+ std::forward_iterator_tag>::value>::type;
+ static constexpr bool DefaultConstructorIsNonTrivial() {
+ return !std::is_trivially_default_constructible<StorageElement>::value;
+ }
+
public:
- // For playing nicely with stl:
- typedef T value_type;
- typedef T* iterator;
- typedef T const* const_iterator;
- typedef T& reference;
- typedef T const& const_reference;
- typedef T* pointer;
- typedef std::ptrdiff_t difference_type;
- typedef size_t size_type;
+ using allocator_type = typename AllocatorTraits::allocator_type;
+ using value_type = typename AllocatorTraits::value_type;
+ using pointer = typename AllocatorTraits::pointer;
+ using const_pointer = typename AllocatorTraits::const_pointer;
+ using reference = value_type&;
+ using const_reference = const value_type&;
+ using size_type = typename AllocatorTraits::size_type;
+ using difference_type = typename AllocatorTraits::difference_type;
+ using iterator = pointer;
+ using const_iterator = const_pointer;
+ using reverse_iterator = std::reverse_iterator<iterator>;
+ using const_reverse_iterator = std::reverse_iterator<const_iterator>;
- // REQUIRES: n >= 0
- // Creates an array object that can store "n" elements.
+ static constexpr size_type inline_elements =
+ (N == kFixedArrayUseDefault ? kInlineBytesDefault / sizeof(value_type)
+ : static_cast<size_type>(N));
+
+ FixedArray(const FixedArray& other,
+ const allocator_type& a = allocator_type())
+ : FixedArray(other.begin(), other.end(), a) {}
+
+ FixedArray(FixedArray&& other, const allocator_type& a = allocator_type())
+ : FixedArray(std::make_move_iterator(other.begin()),
+ std::make_move_iterator(other.end()),
+ a) {}
+
+ // Creates an array object that can store `n` elements.
+ // Note that trivially constructible elements will be uninitialized.
+ explicit FixedArray(size_type n, const allocator_type& a = allocator_type())
+ : storage_(n, a) {
+ if (DefaultConstructorIsNonTrivial()) {
+ ConstructRange(storage_.alloc(), storage_.begin(), storage_.end());
+ }
+ }
+
+ // Creates an array initialized with `n` copies of `val`.
+ FixedArray(size_type n,
+ const value_type& val,
+ const allocator_type& a = allocator_type())
+ : storage_(n, a) {
+ ConstructRange(storage_.alloc(), storage_.begin(), storage_.end(), val);
+ }
+
+ // Creates an array initialized with the size and contents of `init_list`.
+ FixedArray(std::initializer_list<value_type> init_list,
+ const allocator_type& a = allocator_type())
+ : FixedArray(init_list.begin(), init_list.end(), a) {}
+
+ // Creates an array initialized with the elements from the input
+ // range. The array's size will always be `std::distance(first, last)`.
+ // REQUIRES: Iterator must be a forward_iterator or better.
+ template <typename Iterator, EnableIfForwardIterator<Iterator>* = nullptr>
+ FixedArray(Iterator first,
+ Iterator last,
+ const allocator_type& a = allocator_type())
+ : storage_(std::distance(first, last), a) {
+ CopyRange(storage_.alloc(), storage_.begin(), first, last);
+ }
+
+ ~FixedArray() noexcept {
+ for (auto* cur = storage_.begin(); cur != storage_.end(); ++cur) {
+ AllocatorTraits::destroy(storage_.alloc(), cur);
+ }
+ }
+
+ // Assignments are deleted because they break the invariant that the size of a
+ // `FixedArray` never changes.
+ void operator=(FixedArray&&) = delete;
+ void operator=(const FixedArray&) = delete;
+
+ // FixedArray::size()
//
- // FixedArray<T> will not zero-initialize POD (simple) types like int,
- // double, bool, etc.
- // Non-POD types will be default-initialized just like regular vectors or
- // arrays.
- explicit FixedArray(size_type n);
+ // Returns the length of the fixed array.
+ size_type size() const { return storage_.size(); }
- // Releases any resources.
- ~FixedArray();
-
- // Returns the length of the array.
- inline size_type size() const { return size_; }
-
- // Returns the memory size of the array in bytes.
- inline size_t memsize() const { return size_ * sizeof(T); }
-
- // Returns a pointer to the underlying element array.
- inline const T* get() const { return &array_[0].element; }
- inline T* get() { return &array_[0].element; }
-
- // REQUIRES: 0 <= i < size()
- // Returns a reference to the "i"th element.
- inline T& operator[](size_type i) {
- DCHECK_LT(i, size_);
- return array_[i].element;
+ // FixedArray::max_size()
+ //
+ // Returns the largest possible value of `std::distance(begin(), end())` for a
+ // `FixedArray<T>`. This is equivalent to the most possible addressable bytes
+ // over the number of bytes taken by T.
+ constexpr size_type max_size() const {
+ return (std::numeric_limits<difference_type>::max)() / sizeof(value_type);
}
+ // FixedArray::empty()
+ //
+ // Returns whether or not the fixed array is empty.
+ bool empty() const { return size() == 0; }
+
+ // FixedArray::memsize()
+ //
+ // Returns the memory size of the fixed array in bytes.
+ size_t memsize() const { return size() * sizeof(value_type); }
+
+ // FixedArray::data()
+ //
+ // Returns a const T* pointer to elements of the `FixedArray`. This pointer
+ // can be used to access (but not modify) the contained elements.
+ const_pointer data() const { return AsValueType(storage_.begin()); }
+
+ // Overload of FixedArray::data() to return a T* pointer to elements of the
+ // fixed array. This pointer can be used to access and modify the contained
+ // elements.
+ pointer data() { return AsValueType(storage_.begin()); }
+
+ // FixedArray::operator[]
+ //
+ // Returns a reference the ith element of the fixed array.
// REQUIRES: 0 <= i < size()
- // Returns a reference to the "i"th element.
- inline const T& operator[](size_type i) const {
- DCHECK_LT(i, size_);
- return array_[i].element;
+ reference operator[](size_type i) {
+ DCHECK_LT(i, size());
+ return data()[i];
}
- inline iterator begin() { return &array_[0].element; }
- inline iterator end() { return &array_[size_].element; }
+ // Overload of FixedArray::operator()[] to return a const reference to the
+ // ith element of the fixed array.
+ // REQUIRES: 0 <= i < size()
+ const_reference operator[](size_type i) const {
+ DCHECK_LT(i, size());
+ return data()[i];
+ }
- inline const_iterator begin() const { return &array_[0].element; }
- inline const_iterator end() const { return &array_[size_].element; }
+ // FixedArray::front()
+ //
+ // Returns a reference to the first element of the fixed array.
+ reference front() { return *begin(); }
+
+ // Overload of FixedArray::front() to return a reference to the first element
+ // of a fixed array of const values.
+ const_reference front() const { return *begin(); }
+
+ // FixedArray::back()
+ //
+ // Returns a reference to the last element of the fixed array.
+ reference back() { return *(end() - 1); }
+
+ // Overload of FixedArray::back() to return a reference to the last element
+ // of a fixed array of const values.
+ const_reference back() const { return *(end() - 1); }
+
+ // FixedArray::begin()
+ //
+ // Returns an iterator to the beginning of the fixed array.
+ iterator begin() { return data(); }
+
+ // Overload of FixedArray::begin() to return a const iterator to the
+ // beginning of the fixed array.
+ const_iterator begin() const { return data(); }
+
+ // FixedArray::cbegin()
+ //
+ // Returns a const iterator to the beginning of the fixed array.
+ const_iterator cbegin() const { return begin(); }
+
+ // FixedArray::end()
+ //
+ // Returns an iterator to the end of the fixed array.
+ iterator end() { return data() + size(); }
+
+ // Overload of FixedArray::end() to return a const iterator to the end of the
+ // fixed array.
+ const_iterator end() const { return data() + size(); }
+
+ // FixedArray::cend()
+ //
+ // Returns a const iterator to the end of the fixed array.
+ const_iterator cend() const { return end(); }
+
+ // FixedArray::rbegin()
+ //
+ // Returns a reverse iterator from the end of the fixed array.
+ reverse_iterator rbegin() { return reverse_iterator(end()); }
+
+ // Overload of FixedArray::rbegin() to return a const reverse iterator from
+ // the end of the fixed array.
+ const_reverse_iterator rbegin() const {
+ return const_reverse_iterator(end());
+ }
+
+ // FixedArray::crbegin()
+ //
+ // Returns a const reverse iterator from the end of the fixed array.
+ const_reverse_iterator crbegin() const { return rbegin(); }
+
+ // FixedArray::rend()
+ //
+ // Returns a reverse iterator from the beginning of the fixed array.
+ reverse_iterator rend() { return reverse_iterator(begin()); }
+
+ // Overload of FixedArray::rend() for returning a const reverse iterator
+ // from the beginning of the fixed array.
+ const_reverse_iterator rend() const {
+ return const_reverse_iterator(begin());
+ }
+
+ // FixedArray::crend()
+ //
+ // Returns a reverse iterator from the beginning of the fixed array.
+ const_reverse_iterator crend() const { return rend(); }
+
+ // FixedArray::fill()
+ //
+ // Assigns the given `value` to all elements in the fixed array.
+ void fill(const value_type& val) { std::fill(begin(), end(), val); }
+
+ // Relational operators. Equality operators are elementwise using
+ // `operator==`, while order operators order FixedArrays lexicographically.
+ friend bool operator==(const FixedArray& lhs, const FixedArray& rhs) {
+ return std::equal(lhs.begin(), lhs.end(), rhs.begin(), rhs.end());
+ }
+
+ friend bool operator!=(const FixedArray& lhs, const FixedArray& rhs) {
+ return !(lhs == rhs);
+ }
+
+ friend bool operator<(const FixedArray& lhs, const FixedArray& rhs) {
+ return std::lexicographical_compare(
+ lhs.begin(), lhs.end(), rhs.begin(), rhs.end());
+ }
+
+ friend bool operator>(const FixedArray& lhs, const FixedArray& rhs) {
+ return rhs < lhs;
+ }
+
+ friend bool operator<=(const FixedArray& lhs, const FixedArray& rhs) {
+ return !(rhs < lhs);
+ }
+
+ friend bool operator>=(const FixedArray& lhs, const FixedArray& rhs) {
+ return !(lhs < rhs);
+ }
private:
- // Container to hold elements of type T. This is necessary to handle
- // the case where T is a (C-style) array. The size of InnerContainer
- // and T must be the same, otherwise callers' assumptions about use
- // of this code will be broken.
- struct InnerContainer {
- T element;
+ // StorageElement
+ //
+ // For FixedArrays with a C-style-array value_type, StorageElement is a POD
+ // wrapper struct called StorageElementWrapper that holds the value_type
+ // instance inside. This is needed for construction and destruction of the
+ // entire array regardless of how many dimensions it has. For all other cases,
+ // StorageElement is just an alias of value_type.
+ //
+ // Maintainer's Note: The simpler solution would be to simply wrap value_type
+ // in a struct whether it's an array or not. That causes some paranoid
+ // diagnostics to misfire, believing that 'data()' returns a pointer to a
+ // single element, rather than the packed array that it really is.
+ // e.g.:
+ //
+ // FixedArray<char> buf(1);
+ // sprintf(buf.data(), "foo");
+ //
+ // error: call to int __builtin___sprintf_chk(etc...)
+ // will always overflow destination buffer [-Werror]
+ //
+ template <typename OuterT,
+ typename InnerT = typename std::remove_extent<OuterT>::type,
+ size_t InnerN = std::extent<OuterT>::value>
+ struct StorageElementWrapper {
+ InnerT array[InnerN];
};
- // How many elements should we store inline?
- // a. If not specified, use a default of 256 bytes (256 bytes
- // seems small enough to not cause stack overflow or unnecessary
- // stack pollution, while still allowing stack allocation for
- // reasonably long character arrays.
- // b. Never use 0 length arrays (not ISO C++)
- static const size_type S1 = ((inline_elements < 0)
- ? (256/sizeof(T)) : inline_elements);
- static const size_type S2 = (S1 <= 0) ? 1 : S1;
- static const size_type kInlineElements = S2;
+ using StorageElement =
+ typename std::conditional<std::is_array<value_type>::value,
+ StorageElementWrapper<value_type>,
+ value_type>::type;
- size_type const size_;
- InnerContainer* const array_;
+ static pointer AsValueType(pointer ptr) { return ptr; }
+ static pointer AsValueType(StorageElementWrapper<value_type>* ptr) {
+ return std::addressof(ptr->array);
+ }
- // Allocate some space, not an array of elements of type T, so that we can
- // skip calling the T constructors and destructors for space we never use.
- ManualConstructor<InnerContainer> inline_space_[kInlineElements];
+ static_assert(sizeof(StorageElement) == sizeof(value_type), "");
+ static_assert(alignof(StorageElement) == alignof(value_type), "");
+
+ class NonEmptyInlinedStorage {
+ public:
+ StorageElement* data() { return reinterpret_cast<StorageElement*>(buff_); }
+ void AnnotateConstruct(size_type) {}
+ void AnnotateDestruct(size_type) {}
+
+ // #ifdef ADDRESS_SANITIZER
+ // void* RedzoneBegin() { return &redzone_begin_; }
+ // void* RedzoneEnd() { return &redzone_end_ + 1; }
+ // #endif // ADDRESS_SANITIZER
+
+ private:
+ // ADDRESS_SANITIZER_REDZONE(redzone_begin_);
+ alignas(StorageElement) char buff_[sizeof(StorageElement[inline_elements])];
+ // ADDRESS_SANITIZER_REDZONE(redzone_end_);
+ };
+
+ class EmptyInlinedStorage {
+ public:
+ StorageElement* data() { return nullptr; }
+ void AnnotateConstruct(size_type) {}
+ void AnnotateDestruct(size_type) {}
+ };
+
+ using InlinedStorage =
+ typename std::conditional<inline_elements == 0,
+ EmptyInlinedStorage,
+ NonEmptyInlinedStorage>::type;
+
+ // Storage
+ //
+ // An instance of Storage manages the inline and out-of-line memory for
+ // instances of FixedArray. This guarantees that even when construction of
+ // individual elements fails in the FixedArray constructor body, the
+ // destructor for Storage will still be called and out-of-line memory will be
+ // properly deallocated.
+ //
+ class Storage : public InlinedStorage {
+ public:
+ Storage(size_type n, const allocator_type& a)
+ : size_alloc_(n, a), data_(InitializeData()) {}
+
+ ~Storage() noexcept {
+ if (UsingInlinedStorage(size())) {
+ InlinedStorage::AnnotateDestruct(size());
+ } else {
+ AllocatorTraits::deallocate(alloc(), AsValueType(begin()), size());
+ }
+ }
+
+ size_type size() const { return std::get<0>(size_alloc_); }
+ StorageElement* begin() const { return data_; }
+ StorageElement* end() const { return begin() + size(); }
+ allocator_type& alloc() { return std::get<1>(size_alloc_); }
+
+ private:
+ static bool UsingInlinedStorage(size_type n) {
+ return n <= inline_elements;
+ }
+
+ StorageElement* InitializeData() {
+ if (UsingInlinedStorage(size())) {
+ InlinedStorage::AnnotateConstruct(size());
+ return InlinedStorage::data();
+ } else {
+ return reinterpret_cast<StorageElement*>(
+ AllocatorTraits::allocate(alloc(), size()));
+ }
+ }
+
+ // Using std::tuple and not absl::CompressedTuple, as it has a lot of
+ // dependencies to other absl headers.
+ std::tuple<size_type, allocator_type> size_alloc_;
+ StorageElement* data_;
+ };
+
+ Storage storage_;
};
-// Implementation details follow
+template <typename T, size_t N, typename A>
+constexpr size_t FixedArray<T, N, A>::kInlineBytesDefault;
-template <class T, ssize_t S>
-inline FixedArray<T, S>::FixedArray(typename FixedArray<T, S>::size_type n)
- : size_(n),
- array_((n <= kInlineElements
- ? reinterpret_cast<InnerContainer*>(inline_space_)
- : new InnerContainer[n])) {
- // Construct only the elements actually used.
- if (array_ == reinterpret_cast<InnerContainer*>(inline_space_)) {
- for (size_t i = 0; i != size_; ++i) {
- inline_space_[i].Init();
- }
- }
-}
-
-template <class T, ssize_t S>
-inline FixedArray<T, S>::~FixedArray() {
- if (array_ != reinterpret_cast<InnerContainer*>(inline_space_)) {
- delete[] array_;
- } else {
- for (size_t i = 0; i != size_; ++i) {
- inline_space_[i].Destroy();
- }
- }
-}
+template <typename T, size_t N, typename A>
+constexpr typename FixedArray<T, N, A>::size_type
+ FixedArray<T, N, A>::inline_elements;
} // namespace internal
} // namespace ceres
diff --git a/include/ceres/internal/householder_vector.h b/include/ceres/internal/householder_vector.h
new file mode 100644
index 0000000..55f68e5
--- /dev/null
+++ b/include/ceres/internal/householder_vector.h
@@ -0,0 +1,88 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2015 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: vitus@google.com (Michael Vitus)
+
+#ifndef CERES_PUBLIC_INTERNAL_HOUSEHOLDER_VECTOR_H_
+#define CERES_PUBLIC_INTERNAL_HOUSEHOLDER_VECTOR_H_
+
+#include "Eigen/Core"
+#include "glog/logging.h"
+
+namespace ceres {
+namespace internal {
+
+// Algorithm 5.1.1 from 'Matrix Computations' by Golub et al. (Johns Hopkins
+// Studies in Mathematical Sciences) but using the nth element of the input
+// vector as pivot instead of first. This computes the vector v with v(n) = 1
+// and beta such that H = I - beta * v * v^T is orthogonal and
+// H * x = ||x||_2 * e_n.
+//
+// NOTE: Some versions of MSVC have trouble deducing the type of v if
+// you do not specify all the template arguments explicitly.
+template <typename XVectorType, typename Scalar, int N>
+void ComputeHouseholderVector(const XVectorType& x,
+ Eigen::Matrix<Scalar, N, 1>* v,
+ Scalar* beta) {
+ CHECK(beta != nullptr);
+ CHECK(v != nullptr);
+ CHECK_GT(x.rows(), 1);
+ CHECK_EQ(x.rows(), v->rows());
+
+ Scalar sigma = x.head(x.rows() - 1).squaredNorm();
+ *v = x;
+ (*v)(v->rows() - 1) = Scalar(1.0);
+
+ *beta = Scalar(0.0);
+ const Scalar& x_pivot = x(x.rows() - 1);
+
+ if (sigma <= Scalar(std::numeric_limits<double>::epsilon())) {
+ if (x_pivot < Scalar(0.0)) {
+ *beta = Scalar(2.0);
+ }
+ return;
+ }
+
+ const Scalar mu = sqrt(x_pivot * x_pivot + sigma);
+ Scalar v_pivot = Scalar(1.0);
+
+ if (x_pivot <= Scalar(0.0)) {
+ v_pivot = x_pivot - mu;
+ } else {
+ v_pivot = -sigma / (x_pivot + mu);
+ }
+
+ *beta = Scalar(2.0) * v_pivot * v_pivot / (sigma + v_pivot * v_pivot);
+
+ v->head(v->rows() - 1) /= v_pivot;
+}
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_PUBLIC_INTERNAL_HOUSEHOLDER_VECTOR_H_
diff --git a/include/ceres/internal/integer_sequence.h b/include/ceres/internal/integer_sequence.h
deleted file mode 100644
index 3936b92..0000000
--- a/include/ceres/internal/integer_sequence.h
+++ /dev/null
@@ -1,106 +0,0 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2018 Google Inc. All rights reserved.
-// http://ceres-solver.org/
-//
-// 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: jodebo_beck@gmx.de (Johannes Beck)
-//
-// This class mimics std::integer_sequence. That is the reason to follow the
-// naming convention of the stl and not the google one. Once Ceres switches
-// to c++ 14 this class can be removed.
-
-#ifndef CERES_PUBLIC_INTERNAL_INTEGER_SEQUENCE_H_
-#define CERES_PUBLIC_INTERNAL_INTEGER_SEQUENCE_H_
-
-#if __cplusplus >= 201402L
-// We have at least c++ 14 support. Use integer_sequence from the standard.
-// Sometimes the STL implementation uses a compiler intrinsic to generate
-// the sequences which will speed up compilation.
-#include <utility>
-
-namespace ceres {
-namespace internal {
-template <typename T, T... Ns>
-using integer_sequence = std::integer_sequence<T, Ns...>;
-
-template <typename T, T N>
-using make_integer_sequence = std::make_integer_sequence<T, N>;
-
-} // namespace internal
-} // namespace ceres
-#else
-
-namespace ceres {
-namespace internal {
-
-template <typename T, T... Ns>
-struct integer_sequence {
- using value_type = T;
-};
-
-// Implementation of make_integer_sequence.
-//
-// Recursively instantiate make_integer_sequence_impl until Ns
-// contains the sequence 0, 1, ..., Total-1.
-//
-// Example for Total = 4:
-// T CurIdx, Total, Ns...
-// make_integer_sequence_impl<int, 0, 4 >
-// make_integer_sequence_impl<int, 1, 4, 0 >
-// make_integer_sequence_impl<int, 2, 4, 0, 1 >
-// make_integer_sequence_impl<int, 3, 4, 0, 1, 2 >
-// make_integer_sequence_impl<int, 4, 4, 0, 1, 2, 3>
-// ^^^^^^^^^^
-// resulting sequence.
-//
-// The implemented algorithm has linear complexity for simplicity. A O(log(N))
-// implementation can be found e.g. here:
-// https://stackoverflow.com/questions/17424477/implementation-c14-make-integer-sequence
-template <typename T, T CurIdx, T Total, T... Ns>
-struct make_integer_sequence_impl {
- using type = typename make_integer_sequence_impl<T, CurIdx + 1, Total, Ns...,
- CurIdx>::type;
-};
-
-// End of 'recursion' when CurIdx reaches Total. All indices 0, 1, ..., N-1 are
-// contained in Ns. The final integer_sequence is created here.
-template <typename T, T Total, T... Ns>
-struct make_integer_sequence_impl<T, Total, Total, Ns...> {
- using type = integer_sequence<T, Ns...>;
-};
-
-// A helper alias template to simplify creation of integer_sequence with 0, 1,
-// ..., N-1 as Ns:
-template <typename T, T N>
-using make_integer_sequence =
- typename make_integer_sequence_impl<T, 0, N>::type;
-
-} // namespace internal
-} // namespace ceres
-
-#endif
-
-#endif // CERES_PUBLIC_INTERNAL_INTEGER_SEQUENCE_H_
diff --git a/include/ceres/internal/integer_sequence_algorithm.h b/include/ceres/internal/integer_sequence_algorithm.h
index bdeab93..8c0f3bc 100644
--- a/include/ceres/internal/integer_sequence_algorithm.h
+++ b/include/ceres/internal/integer_sequence_algorithm.h
@@ -35,7 +35,7 @@
#ifndef CERES_PUBLIC_INTERNAL_INTEGER_SEQUENCE_ALGORITHM_H_
#define CERES_PUBLIC_INTERNAL_INTEGER_SEQUENCE_ALGORITHM_H_
-#include "integer_sequence.h"
+#include <utility>
namespace ceres {
namespace internal {
@@ -61,33 +61,34 @@
// Strip of and sum the first number.
template <typename T, T N, T... Ns>
-struct SumImpl<integer_sequence<T, N, Ns...>> {
- static constexpr T Value = N + SumImpl<integer_sequence<T, Ns...>>::Value;
+struct SumImpl<std::integer_sequence<T, N, Ns...>> {
+ static constexpr T Value =
+ N + SumImpl<std::integer_sequence<T, Ns...>>::Value;
};
// Strip of and sum the first two numbers.
template <typename T, T N1, T N2, T... Ns>
-struct SumImpl<integer_sequence<T, N1, N2, Ns...>> {
+struct SumImpl<std::integer_sequence<T, N1, N2, Ns...>> {
static constexpr T Value =
- N1 + N2 + SumImpl<integer_sequence<T, Ns...>>::Value;
+ N1 + N2 + SumImpl<std::integer_sequence<T, Ns...>>::Value;
};
// Strip of and sum the first four numbers.
template <typename T, T N1, T N2, T N3, T N4, T... Ns>
-struct SumImpl<integer_sequence<T, N1, N2, N3, N4, Ns...>> {
+struct SumImpl<std::integer_sequence<T, N1, N2, N3, N4, Ns...>> {
static constexpr T Value =
- N1 + N2 + N3 + N4 + SumImpl<integer_sequence<T, Ns...>>::Value;
+ N1 + N2 + N3 + N4 + SumImpl<std::integer_sequence<T, Ns...>>::Value;
};
// Only one number is left. 'Value' is just that number ('recursion' ends).
template <typename T, T N>
-struct SumImpl<integer_sequence<T, N>> {
+struct SumImpl<std::integer_sequence<T, N>> {
static constexpr T Value = N;
};
// No number is left. 'Value' is the identity element (for sum this is zero).
template <typename T>
-struct SumImpl<integer_sequence<T>> {
+struct SumImpl<std::integer_sequence<T>> {
static constexpr T Value = T(0);
};
@@ -129,16 +130,20 @@
struct ExclusiveScanImpl;
template <typename T, T Sum, T N, T... Ns, T... Rs>
-struct ExclusiveScanImpl<T, Sum, integer_sequence<T, N, Ns...>,
- integer_sequence<T, Rs...>> {
+struct ExclusiveScanImpl<T,
+ Sum,
+ std::integer_sequence<T, N, Ns...>,
+ std::integer_sequence<T, Rs...>> {
using Type =
- typename ExclusiveScanImpl<T, Sum + N, integer_sequence<T, Ns...>,
- integer_sequence<T, Rs..., Sum>>::Type;
+ typename ExclusiveScanImpl<T,
+ Sum + N,
+ std::integer_sequence<T, Ns...>,
+ std::integer_sequence<T, Rs..., Sum>>::Type;
};
// End of 'recursion'. The resulting type is SeqOut.
template <typename T, T Sum, typename SeqOut>
-struct ExclusiveScanImpl<T, Sum, integer_sequence<T>, SeqOut> {
+struct ExclusiveScanImpl<T, Sum, std::integer_sequence<T>, SeqOut> {
using Type = SeqOut;
};
@@ -152,7 +157,7 @@
public:
using Type =
- typename ExclusiveScanImpl<T, T(0), Seq, integer_sequence<T>>::Type;
+ typename ExclusiveScanImpl<T, T(0), Seq, std::integer_sequence<T>>::Type;
};
// Helper to use exclusive scan without typename.
diff --git a/include/ceres/internal/line_parameterization.h b/include/ceres/internal/line_parameterization.h
new file mode 100644
index 0000000..eda3901
--- /dev/null
+++ b/include/ceres/internal/line_parameterization.h
@@ -0,0 +1,183 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2020 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// 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: jodebo_beck@gmx.de (Johannes Beck)
+//
+
+#ifndef CERES_PUBLIC_INTERNAL_LINE_PARAMETERIZATION_H_
+#define CERES_PUBLIC_INTERNAL_LINE_PARAMETERIZATION_H_
+
+#include "householder_vector.h"
+
+namespace ceres {
+
+template <int AmbientSpaceDimension>
+bool LineParameterization<AmbientSpaceDimension>::Plus(
+ const double* x_ptr,
+ const double* delta_ptr,
+ double* x_plus_delta_ptr) const {
+ // We seek a box plus operator of the form
+ //
+ // [o*, d*] = Plus([o, d], [delta_o, delta_d])
+ //
+ // where o is the origin point, d is the direction vector, delta_o is
+ // the delta of the origin point and delta_d the delta of the direction and
+ // o* and d* is the updated origin point and direction.
+ //
+ // We separate the Plus operator into the origin point and directional part
+ // d* = Plus_d(d, delta_d)
+ // o* = Plus_o(o, d, delta_o)
+ //
+ // The direction update function Plus_d is the same as for the homogeneous
+ // vector parameterization:
+ //
+ // d* = H_{v(d)} [0.5 sinc(0.5 |delta_d|) delta_d, cos(0.5 |delta_d|)]^T
+ //
+ // where H is the householder matrix
+ // H_{v} = I - (2 / |v|^2) v v^T
+ // and
+ // v(d) = d - sign(d_n) |d| e_n.
+ //
+ // The origin point update function Plus_o is defined as
+ //
+ // o* = o + H_{v(d)} [0.5 delta_o, 0]^T.
+
+ static constexpr int kDim = AmbientSpaceDimension;
+ using AmbientVector = Eigen::Matrix<double, kDim, 1>;
+ using AmbientVectorRef = Eigen::Map<Eigen::Matrix<double, kDim, 1>>;
+ using ConstAmbientVectorRef =
+ Eigen::Map<const Eigen::Matrix<double, kDim, 1>>;
+ using ConstTangentVectorRef =
+ Eigen::Map<const Eigen::Matrix<double, kDim - 1, 1>>;
+
+ ConstAmbientVectorRef o(x_ptr);
+ ConstAmbientVectorRef d(x_ptr + kDim);
+
+ ConstTangentVectorRef delta_o(delta_ptr);
+ ConstTangentVectorRef delta_d(delta_ptr + kDim - 1);
+ AmbientVectorRef o_plus_delta(x_plus_delta_ptr);
+ AmbientVectorRef d_plus_delta(x_plus_delta_ptr + kDim);
+
+ const double norm_delta_d = delta_d.norm();
+
+ o_plus_delta = o;
+
+ // Shortcut for zero delta direction.
+ if (norm_delta_d == 0.0) {
+ d_plus_delta = d;
+
+ if (delta_o.isZero(0.0)) {
+ return true;
+ }
+ }
+
+ // Calculate the householder transformation which is needed for f_d and f_o.
+ AmbientVector v;
+ double beta;
+
+ // NOTE: The explicit template arguments are needed here because
+ // ComputeHouseholderVector is templated and some versions of MSVC
+ // have trouble deducing the type of v automatically.
+ internal::ComputeHouseholderVector<ConstAmbientVectorRef, double, kDim>(
+ d, &v, &beta);
+
+ if (norm_delta_d != 0.0) {
+ // Map the delta from the minimum representation to the over parameterized
+ // homogeneous vector. See section A6.9.2 on page 624 of Hartley & Zisserman
+ // (2nd Edition) for a detailed description. Note there is a typo on Page
+ // 625, line 4 so check the book errata.
+ const double norm_delta_div_2 = 0.5 * norm_delta_d;
+ const double sin_delta_by_delta =
+ std::sin(norm_delta_div_2) / norm_delta_div_2;
+
+ // Apply the delta update to remain on the unit sphere. See section A6.9.3
+ // on page 625 of Hartley & Zisserman (2nd Edition) for a detailed
+ // description.
+ AmbientVector y;
+ y.template head<kDim - 1>() = 0.5 * sin_delta_by_delta * delta_d;
+ y[kDim - 1] = std::cos(norm_delta_div_2);
+
+ d_plus_delta = d.norm() * (y - v * (beta * (v.transpose() * y)));
+ }
+
+ // The null space is in the direction of the line, so the tangent space is
+ // perpendicular to the line direction. This is achieved by using the
+ // householder matrix of the direction and allow only movements
+ // perpendicular to e_n.
+ //
+ // The factor of 0.5 is used to be consistent with the line direction
+ // update.
+ AmbientVector y;
+ y << 0.5 * delta_o, 0;
+ o_plus_delta += y - v * (beta * (v.transpose() * y));
+
+ return true;
+}
+
+template <int AmbientSpaceDimension>
+bool LineParameterization<AmbientSpaceDimension>::ComputeJacobian(
+ const double* x_ptr, double* jacobian_ptr) const {
+ static constexpr int kDim = AmbientSpaceDimension;
+ using AmbientVector = Eigen::Matrix<double, kDim, 1>;
+ using ConstAmbientVectorRef =
+ Eigen::Map<const Eigen::Matrix<double, kDim, 1>>;
+ using MatrixRef = Eigen::Map<
+ Eigen::Matrix<double, 2 * kDim, 2 * (kDim - 1), Eigen::RowMajor>>;
+
+ ConstAmbientVectorRef d(x_ptr + kDim);
+ MatrixRef jacobian(jacobian_ptr);
+
+ // Clear the Jacobian as only half of the matrix is not zero.
+ jacobian.setZero();
+
+ AmbientVector v;
+ double beta;
+
+ // NOTE: The explicit template arguments are needed here because
+ // ComputeHouseholderVector is templated and some versions of MSVC
+ // have trouble deducing the type of v automatically.
+ internal::ComputeHouseholderVector<ConstAmbientVectorRef, double, kDim>(
+ d, &v, &beta);
+
+ // The Jacobian is equal to J = 0.5 * H.leftCols(kDim - 1) where H is
+ // the Householder matrix (H = I - beta * v * v') for the origin point. For
+ // the line direction part the Jacobian is scaled by the norm of the
+ // direction.
+ for (int i = 0; i < kDim - 1; ++i) {
+ jacobian.block(0, i, kDim, 1) = -0.5 * beta * v(i) * v;
+ jacobian.col(i)(i) += 0.5;
+ }
+
+ jacobian.template block<kDim, kDim - 1>(kDim, kDim - 1) =
+ jacobian.template block<kDim, kDim - 1>(0, 0) * d.norm();
+ return true;
+}
+
+} // namespace ceres
+
+#endif // CERES_PUBLIC_INTERNAL_LINE_PARAMETERIZATION_H_
diff --git a/include/ceres/internal/manual_constructor.h b/include/ceres/internal/manual_constructor.h
deleted file mode 100644
index 60f147a..0000000
--- a/include/ceres/internal/manual_constructor.h
+++ /dev/null
@@ -1,152 +0,0 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
-// http://ceres-solver.org/
-//
-// 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: kenton@google.com (Kenton Varda)
-//
-// ManualConstructor statically-allocates space in which to store some
-// object, but does not initialize it. You can then call the constructor
-// and destructor for the object yourself as you see fit. This is useful
-// for memory management optimizations, where you want to initialize and
-// destroy an object multiple times but only allocate it once.
-//
-// (When I say ManualConstructor statically allocates space, I mean that
-// the ManualConstructor object itself is forced to be the right size.)
-
-#ifndef CERES_PUBLIC_INTERNAL_MANUAL_CONSTRUCTOR_H_
-#define CERES_PUBLIC_INTERNAL_MANUAL_CONSTRUCTOR_H_
-
-#include <new>
-#include <utility>
-
-namespace ceres {
-namespace internal {
-
-// ------- Define CERES_ALIGNED_CHAR_ARRAY --------------------------------
-
-// Platform independent macros to get aligned memory allocations.
-// For example
-//
-// MyFoo my_foo CERES_ALIGN_ATTRIBUTE(16);
-//
-// Gives us an instance of MyFoo which is aligned at a 16 byte
-// boundary.
-#if defined(_MSC_VER)
-#define CERES_ALIGN_ATTRIBUTE(n) __declspec(align(n))
-#define CERES_ALIGN_OF(T) __alignof(T)
-#elif defined(__GNUC__)
-#define CERES_ALIGN_ATTRIBUTE(n) __attribute__((aligned(n)))
-#define CERES_ALIGN_OF(T) __alignof(T)
-#endif
-
-#ifndef CERES_ALIGNED_CHAR_ARRAY
-
-// Because MSVC and older GCCs require that the argument to their alignment
-// construct to be a literal constant integer, we use a template instantiated
-// at all the possible powers of two.
-template<int alignment, int size> struct AlignType { };
-template<int size> struct AlignType<0, size> { typedef char result[size]; };
-
-#if !defined(CERES_ALIGN_ATTRIBUTE)
-#define CERES_ALIGNED_CHAR_ARRAY you_must_define_CERES_ALIGNED_CHAR_ARRAY_for_your_compiler
-#else // !defined(CERES_ALIGN_ATTRIBUTE)
-
-#define CERES_ALIGN_TYPE_TEMPLATE(X) \
- template<int size> struct AlignType<X, size> { \
- typedef CERES_ALIGN_ATTRIBUTE(X) char result[size]; \
- }
-
-CERES_ALIGN_TYPE_TEMPLATE(1);
-CERES_ALIGN_TYPE_TEMPLATE(2);
-CERES_ALIGN_TYPE_TEMPLATE(4);
-CERES_ALIGN_TYPE_TEMPLATE(8);
-CERES_ALIGN_TYPE_TEMPLATE(16);
-CERES_ALIGN_TYPE_TEMPLATE(32);
-CERES_ALIGN_TYPE_TEMPLATE(64);
-CERES_ALIGN_TYPE_TEMPLATE(128);
-CERES_ALIGN_TYPE_TEMPLATE(256);
-CERES_ALIGN_TYPE_TEMPLATE(512);
-CERES_ALIGN_TYPE_TEMPLATE(1024);
-CERES_ALIGN_TYPE_TEMPLATE(2048);
-CERES_ALIGN_TYPE_TEMPLATE(4096);
-CERES_ALIGN_TYPE_TEMPLATE(8192);
-// Any larger and MSVC++ will complain.
-
-#undef CERES_ALIGN_TYPE_TEMPLATE
-
-#define CERES_ALIGNED_CHAR_ARRAY(T, Size) \
- typename AlignType<CERES_ALIGN_OF(T), sizeof(T) * Size>::result
-
-#endif // !defined(CERES_ALIGN_ATTRIBUTE)
-
-#endif // CERES_ALIGNED_CHAR_ARRAY
-
-template <typename Type>
-class ManualConstructor {
- public:
- // No constructor or destructor because one of the most useful uses of
- // this class is as part of a union, and members of a union cannot have
- // constructors or destructors. And, anyway, the whole point of this
- // class is to bypass these.
-
- inline Type* get() {
- return reinterpret_cast<Type*>(space_);
- }
- inline const Type* get() const {
- return reinterpret_cast<const Type*>(space_);
- }
-
- inline Type* operator->() { return get(); }
- inline const Type* operator->() const { return get(); }
-
- inline Type& operator*() { return *get(); }
- inline const Type& operator*() const { return *get(); }
-
- // This is needed to get around the strict aliasing warning GCC generates.
- inline void* space() {
- return reinterpret_cast<void*>(space_);
- }
-
- template <typename... Ts>
- inline void Init(Ts&&... ps) {
- new(space()) Type(std::forward<Ts>(ps)...);
- }
-
- inline void Destroy() {
- get()->~Type();
- }
-
- private:
- CERES_ALIGNED_CHAR_ARRAY(Type, 1) space_;
-};
-
-#undef CERES_ALIGNED_CHAR_ARRAY
-
-} // namespace internal
-} // namespace ceres
-
-#endif // CERES_PUBLIC_INTERNAL_MANUAL_CONSTRUCTOR_H_
diff --git a/include/ceres/internal/memory.h b/include/ceres/internal/memory.h
new file mode 100644
index 0000000..45c5b67
--- /dev/null
+++ b/include/ceres/internal/memory.h
@@ -0,0 +1,90 @@
+// Copyright 2017 The Abseil Authors.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+//
+// -----------------------------------------------------------------------------
+// File: memory.h
+// -----------------------------------------------------------------------------
+//
+// This header file contains utility functions for managing the creation and
+// conversion of smart pointers. This file is an extension to the C++
+// standard <memory> library header file.
+
+#ifndef CERES_PUBLIC_INTERNAL_MEMORY_H_
+#define CERES_PUBLIC_INTERNAL_MEMORY_H_
+
+#include <memory>
+
+#ifdef CERES_HAVE_EXCEPTIONS
+#define CERES_INTERNAL_TRY try
+#define CERES_INTERNAL_CATCH_ANY catch (...)
+#define CERES_INTERNAL_RETHROW \
+ do { \
+ throw; \
+ } while (false)
+#else // CERES_HAVE_EXCEPTIONS
+#define CERES_INTERNAL_TRY if (true)
+#define CERES_INTERNAL_CATCH_ANY else if (false)
+#define CERES_INTERNAL_RETHROW \
+ do { \
+ } while (false)
+#endif // CERES_HAVE_EXCEPTIONS
+
+namespace ceres {
+namespace internal {
+
+template <typename Allocator, typename Iterator, typename... Args>
+void ConstructRange(Allocator& alloc,
+ Iterator first,
+ Iterator last,
+ const Args&... args) {
+ for (Iterator cur = first; cur != last; ++cur) {
+ CERES_INTERNAL_TRY {
+ std::allocator_traits<Allocator>::construct(
+ alloc, std::addressof(*cur), args...);
+ }
+ CERES_INTERNAL_CATCH_ANY {
+ while (cur != first) {
+ --cur;
+ std::allocator_traits<Allocator>::destroy(alloc, std::addressof(*cur));
+ }
+ CERES_INTERNAL_RETHROW;
+ }
+ }
+}
+
+template <typename Allocator, typename Iterator, typename InputIterator>
+void CopyRange(Allocator& alloc,
+ Iterator destination,
+ InputIterator first,
+ InputIterator last) {
+ for (Iterator cur = destination; first != last;
+ static_cast<void>(++cur), static_cast<void>(++first)) {
+ CERES_INTERNAL_TRY {
+ std::allocator_traits<Allocator>::construct(
+ alloc, std::addressof(*cur), *first);
+ }
+ CERES_INTERNAL_CATCH_ANY {
+ while (cur != destination) {
+ --cur;
+ std::allocator_traits<Allocator>::destroy(alloc, std::addressof(*cur));
+ }
+ CERES_INTERNAL_RETHROW;
+ }
+ }
+}
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_PUBLIC_INTERNAL_MEMORY_H_
diff --git a/include/ceres/internal/numeric_diff.h b/include/ceres/internal/numeric_diff.h
index 8851bc7..ff7a2c3 100644
--- a/include/ceres/internal/numeric_diff.h
+++ b/include/ceres/internal/numeric_diff.h
@@ -36,6 +36,7 @@
#define CERES_PUBLIC_INTERNAL_NUMERIC_DIFF_H_
#include <cstring>
+#include <utility>
#include "Eigen/Dense"
#include "Eigen/StdVector"
@@ -46,15 +47,17 @@
#include "ceres/types.h"
#include "glog/logging.h"
-
namespace ceres {
namespace internal {
// 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, NumericDiffMethodType kMethod,
- int kNumResiduals, typename ParameterDims, int kParameterBlock,
+template <typename CostFunctor,
+ NumericDiffMethodType kMethod,
+ int kNumResiduals,
+ typename ParameterDims,
+ int kParameterBlock,
int kParameterBlockSize>
struct NumericDiff {
// Mutates parameters but must restore them before return.
@@ -65,23 +68,23 @@
int num_residuals,
int parameter_block_index,
int parameter_block_size,
- double **parameters,
- double *jacobian) {
+ double** parameters,
+ double* jacobian) {
+ using Eigen::ColMajor;
using Eigen::Map;
using Eigen::Matrix;
using Eigen::RowMajor;
- using Eigen::ColMajor;
DCHECK(jacobian);
const int num_residuals_internal =
(kNumResiduals != ceres::DYNAMIC ? kNumResiduals : num_residuals);
const int parameter_block_index_internal =
- (kParameterBlock != ceres::DYNAMIC ? kParameterBlock :
- parameter_block_index);
+ (kParameterBlock != ceres::DYNAMIC ? kParameterBlock
+ : parameter_block_index);
const int parameter_block_size_internal =
- (kParameterBlockSize != ceres::DYNAMIC ? kParameterBlockSize :
- parameter_block_size);
+ (kParameterBlockSize != ceres::DYNAMIC ? kParameterBlockSize
+ : parameter_block_size);
typedef Matrix<double, kNumResiduals, 1> ResidualVector;
typedef Matrix<double, kParameterBlockSize, 1> ParameterVector;
@@ -96,17 +99,17 @@
(kParameterBlockSize == 1) ? ColMajor : RowMajor>
JacobianMatrix;
- Map<JacobianMatrix> parameter_jacobian(jacobian,
- num_residuals_internal,
- parameter_block_size_internal);
+ Map<JacobianMatrix> parameter_jacobian(
+ jacobian, num_residuals_internal, parameter_block_size_internal);
Map<ParameterVector> x_plus_delta(
parameters[parameter_block_index_internal],
parameter_block_size_internal);
ParameterVector x(x_plus_delta);
- ParameterVector step_size = x.array().abs() *
- ((kMethod == RIDDERS) ? options.ridders_relative_initial_step_size :
- options.relative_step_size);
+ ParameterVector step_size =
+ x.array().abs() * ((kMethod == RIDDERS)
+ ? options.ridders_relative_initial_step_size
+ : options.relative_step_size);
// It is not a good idea to make the step size arbitrarily
// small. This will lead to problems with round off and numerical
@@ -117,22 +120,24 @@
// For Ridders' method, the initial step size is required to be large,
// thus ridders_relative_initial_step_size is used.
if (kMethod == RIDDERS) {
- min_step_size = std::max(min_step_size,
- options.ridders_relative_initial_step_size);
+ min_step_size =
+ std::max(min_step_size, options.ridders_relative_initial_step_size);
}
// For each parameter in the parameter block, use finite differences to
// compute the derivative for that parameter.
FixedArray<double> temp_residual_array(num_residuals_internal);
FixedArray<double> residual_array(num_residuals_internal);
- Map<ResidualVector> residuals(residual_array.get(),
+ Map<ResidualVector> residuals(residual_array.data(),
num_residuals_internal);
for (int j = 0; j < parameter_block_size_internal; ++j) {
const double delta = std::max(min_step_size, step_size(j));
if (kMethod == RIDDERS) {
- if (!EvaluateRiddersJacobianColumn(functor, j, delta,
+ if (!EvaluateRiddersJacobianColumn(functor,
+ j,
+ delta,
options,
num_residuals_internal,
parameter_block_size_internal,
@@ -140,20 +145,22 @@
residuals_at_eval_point,
parameters,
x_plus_delta.data(),
- temp_residual_array.get(),
- residual_array.get())) {
+ temp_residual_array.data(),
+ residual_array.data())) {
return false;
}
} else {
- if (!EvaluateJacobianColumn(functor, j, delta,
+ if (!EvaluateJacobianColumn(functor,
+ j,
+ delta,
num_residuals_internal,
parameter_block_size_internal,
x.data(),
residuals_at_eval_point,
parameters,
x_plus_delta.data(),
- temp_residual_array.get(),
- residual_array.get())) {
+ temp_residual_array.data(),
+ residual_array.data())) {
return false;
}
}
@@ -181,8 +188,7 @@
typedef Matrix<double, kParameterBlockSize, 1> ParameterVector;
Map<const ParameterVector> x(x_ptr, parameter_block_size);
- Map<ParameterVector> x_plus_delta(x_plus_delta_ptr,
- parameter_block_size);
+ Map<ParameterVector> x_plus_delta(x_plus_delta_ptr, parameter_block_size);
Map<ResidualVector> residuals(residuals_ptr, num_residuals);
Map<ResidualVector> temp_residuals(temp_residuals_ptr, num_residuals);
@@ -190,9 +196,8 @@
// Mutate 1 element at a time and then restore.
x_plus_delta(parameter_index) = x(parameter_index) + delta;
- if (!VariadicEvaluate<ParameterDims>(*functor,
- parameters,
- residuals.data())) {
+ if (!VariadicEvaluate<ParameterDims>(
+ *functor, parameters, residuals.data())) {
return false;
}
@@ -205,9 +210,8 @@
// Compute the function on the other side of x(parameter_index).
x_plus_delta(parameter_index) = x(parameter_index) - delta;
- if (!VariadicEvaluate<ParameterDims>(*functor,
- parameters,
- temp_residuals.data())) {
+ if (!VariadicEvaluate<ParameterDims>(
+ *functor, parameters, temp_residuals.data())) {
return false;
}
@@ -216,8 +220,7 @@
} else {
// Forward difference only; reuse existing residuals evaluation.
residuals -=
- Map<const ResidualVector>(residuals_at_eval_point,
- num_residuals);
+ Map<const ResidualVector>(residuals_at_eval_point, num_residuals);
}
// Restore x_plus_delta.
@@ -253,17 +256,17 @@
double* x_plus_delta_ptr,
double* temp_residuals_ptr,
double* residuals_ptr) {
+ using Eigen::aligned_allocator;
using Eigen::Map;
using Eigen::Matrix;
- using Eigen::aligned_allocator;
typedef Matrix<double, kNumResiduals, 1> ResidualVector;
- typedef Matrix<double, kNumResiduals, Eigen::Dynamic> ResidualCandidateMatrix;
+ typedef Matrix<double, kNumResiduals, Eigen::Dynamic>
+ ResidualCandidateMatrix;
typedef Matrix<double, kParameterBlockSize, 1> ParameterVector;
Map<const ParameterVector> x(x_ptr, parameter_block_size);
- Map<ParameterVector> x_plus_delta(x_plus_delta_ptr,
- parameter_block_size);
+ Map<ParameterVector> x_plus_delta(x_plus_delta_ptr, parameter_block_size);
Map<ResidualVector> residuals(residuals_ptr, num_residuals);
Map<ResidualVector> temp_residuals(temp_residuals_ptr, num_residuals);
@@ -274,18 +277,16 @@
// As the derivative is estimated, the step size decreases.
// By default, the step sizes are chosen so that the middle column
// of the Romberg tableau uses the input delta.
- double current_step_size = delta *
- pow(options.ridders_step_shrink_factor,
- options.max_num_ridders_extrapolations / 2);
+ double current_step_size =
+ delta * pow(options.ridders_step_shrink_factor,
+ options.max_num_ridders_extrapolations / 2);
// Double-buffering temporary differential candidate vectors
// from previous step size.
ResidualCandidateMatrix stepsize_candidates_a(
- num_residuals,
- options.max_num_ridders_extrapolations);
+ num_residuals, options.max_num_ridders_extrapolations);
ResidualCandidateMatrix stepsize_candidates_b(
- num_residuals,
- options.max_num_ridders_extrapolations);
+ num_residuals, options.max_num_ridders_extrapolations);
ResidualCandidateMatrix* current_candidates = &stepsize_candidates_a;
ResidualCandidateMatrix* previous_candidates = &stepsize_candidates_b;
@@ -303,7 +304,9 @@
// 3. Extrapolation becomes numerically unstable.
for (int i = 0; i < options.max_num_ridders_extrapolations; ++i) {
// Compute the numerical derivative at this step size.
- if (!EvaluateJacobianColumn(functor, parameter_index, current_step_size,
+ if (!EvaluateJacobianColumn(functor,
+ parameter_index,
+ current_step_size,
num_residuals,
parameter_block_size,
x.data(),
@@ -326,23 +329,24 @@
// Extrapolation factor for Richardson acceleration method (see below).
double richardson_factor = options.ridders_step_shrink_factor *
- options.ridders_step_shrink_factor;
+ options.ridders_step_shrink_factor;
for (int k = 1; k <= i; ++k) {
// Extrapolate the various orders of finite differences using
// the Richardson acceleration method.
current_candidates->col(k) =
(richardson_factor * current_candidates->col(k - 1) -
- previous_candidates->col(k - 1)) / (richardson_factor - 1.0);
+ previous_candidates->col(k - 1)) /
+ (richardson_factor - 1.0);
richardson_factor *= options.ridders_step_shrink_factor *
- options.ridders_step_shrink_factor;
+ options.ridders_step_shrink_factor;
// Compute the difference between the previous value and the current.
double candidate_error = std::max(
- (current_candidates->col(k) -
- current_candidates->col(k - 1)).norm(),
- (current_candidates->col(k) -
- previous_candidates->col(k - 1)).norm());
+ (current_candidates->col(k) - current_candidates->col(k - 1))
+ .norm(),
+ (current_candidates->col(k) - previous_candidates->col(k - 1))
+ .norm());
// If the error has decreased, update results.
if (candidate_error <= norm_error) {
@@ -364,8 +368,9 @@
// Check to see if the current gradient estimate is numerically unstable.
// If so, bail out and return the last stable result.
if (i > 0) {
- double tableau_error = (current_candidates->col(i) -
- previous_candidates->col(i - 1)).norm();
+ double tableau_error =
+ (current_candidates->col(i) - previous_candidates->col(i - 1))
+ .norm();
// Compare current error to the chosen candidate's error.
if (tableau_error >= 2 * norm_error) {
@@ -437,7 +442,7 @@
template <typename ParameterDims, int N, int... Ns, int ParameterIdx>
struct EvaluateJacobianForParameterBlocks<ParameterDims,
- integer_sequence<int, N, Ns...>,
+ std::integer_sequence<int, N, Ns...>,
ParameterIdx> {
template <NumericDiffMethodType method,
int kNumResiduals,
@@ -468,7 +473,7 @@
}
return EvaluateJacobianForParameterBlocks<ParameterDims,
- integer_sequence<int, Ns...>,
+ std::integer_sequence<int, Ns...>,
ParameterIdx + 1>::
template Apply<method, kNumResiduals>(functor,
residuals_at_eval_point,
@@ -481,14 +486,18 @@
// End of 'recursion'. Nothing more to do.
template <typename ParameterDims, int ParameterIdx>
-struct EvaluateJacobianForParameterBlocks<ParameterDims, integer_sequence<int>,
+struct EvaluateJacobianForParameterBlocks<ParameterDims,
+ std::integer_sequence<int>,
ParameterIdx> {
- template <NumericDiffMethodType method, int kNumResiduals,
+ template <NumericDiffMethodType method,
+ int kNumResiduals,
typename CostFunctor>
static bool Apply(const CostFunctor* /* NOT USED*/,
const double* /* NOT USED*/,
- const NumericDiffOptions& /* NOT USED*/, int /* NOT USED*/,
- double** /* NOT USED*/, double** /* NOT USED*/) {
+ const NumericDiffOptions& /* NOT USED*/,
+ int /* NOT USED*/,
+ double** /* NOT USED*/,
+ double** /* NOT USED*/) {
return true;
}
};
diff --git a/include/ceres/internal/parameter_dims.h b/include/ceres/internal/parameter_dims.h
index 2272e8d..2402106 100644
--- a/include/ceres/internal/parameter_dims.h
+++ b/include/ceres/internal/parameter_dims.h
@@ -32,8 +32,8 @@
#define CERES_PUBLIC_INTERNAL_PARAMETER_DIMS_H_
#include <array>
+#include <utility>
-#include "ceres/internal/integer_sequence.h"
#include "ceres/internal/integer_sequence_algorithm.h"
namespace ceres {
@@ -41,16 +41,16 @@
// Checks, whether the given parameter block sizes are valid. Valid means every
// dimension is bigger than zero.
-constexpr bool IsValidParameterDimensionSequence(integer_sequence<int>) {
+constexpr bool IsValidParameterDimensionSequence(std::integer_sequence<int>) {
return true;
}
template <int N, int... Ts>
constexpr bool IsValidParameterDimensionSequence(
- integer_sequence<int, N, Ts...>) {
+ std::integer_sequence<int, N, Ts...>) {
return (N <= 0) ? false
: IsValidParameterDimensionSequence(
- integer_sequence<int, Ts...>());
+ std::integer_sequence<int, Ts...>());
}
// Helper class that represents the parameter dimensions. The parameter
@@ -66,7 +66,7 @@
template <bool IsDynamic, int... Ns>
class ParameterDims {
public:
- using Parameters = integer_sequence<int, Ns...>;
+ using Parameters = std::integer_sequence<int, Ns...>;
// The parameter dimensions are only valid if all parameter block dimensions
// are greater than zero.
@@ -82,7 +82,7 @@
"At least one parameter block must be specified.");
static constexpr int kNumParameters =
- Sum<integer_sequence<int, Ns...>>::Value;
+ Sum<std::integer_sequence<int, Ns...>>::Value;
static constexpr int GetDim(int dim) { return params_[dim]; }
@@ -98,7 +98,7 @@
private:
template <typename T, int... Indices>
static inline std::array<T*, kNumParameterBlocks> GetUnpackedParameters(
- T* ptr, integer_sequence<int, Indices...>) {
+ T* ptr, std::integer_sequence<int, Indices...>) {
return std::array<T*, kNumParameterBlocks>{{ptr + Indices...}};
}
diff --git a/include/ceres/internal/port.h b/include/ceres/internal/port.h
index 07bb505..040a1ef 100644
--- a/include/ceres/internal/port.h
+++ b/include/ceres/internal/port.h
@@ -35,56 +35,76 @@
#include "ceres/internal/config.h"
#if defined(CERES_USE_OPENMP)
-# if defined(CERES_USE_CXX11_THREADS) || defined(CERES_NO_THREADS)
-# error CERES_USE_OPENMP is mutually exclusive to CERES_USE_CXX11_THREADS and CERES_NO_THREADS
-# endif
-#elif defined(CERES_USE_CXX11_THREADS)
-# if defined(CERES_USE_OPENMP) || defined(CERES_NO_THREADS)
-# error CERES_USE_CXX11_THREADS is mutually exclusive to CERES_USE_OPENMP, CERES_USE_CXX11_THREADS and CERES_NO_THREADS
-# endif
+#if defined(CERES_USE_CXX_THREADS) || defined(CERES_NO_THREADS)
+#error CERES_USE_OPENMP is mutually exclusive to CERES_USE_CXX_THREADS and CERES_NO_THREADS
+#endif
+#elif defined(CERES_USE_CXX_THREADS)
+#if defined(CERES_USE_OPENMP) || defined(CERES_NO_THREADS)
+#error CERES_USE_CXX_THREADS is mutually exclusive to CERES_USE_OPENMP, CERES_USE_CXX_THREADS and CERES_NO_THREADS
+#endif
#elif defined(CERES_NO_THREADS)
-# if defined(CERES_USE_OPENMP) || defined(CERES_USE_CXX11_THREADS)
-# error CERES_NO_THREADS is mutually exclusive to CERES_USE_OPENMP and CERES_USE_CXX11_THREADS
-# endif
+#if defined(CERES_USE_OPENMP) || defined(CERES_USE_CXX_THREADS)
+#error CERES_NO_THREADS is mutually exclusive to CERES_USE_OPENMP and CERES_USE_CXX_THREADS
+#endif
#else
-# error One of CERES_USE_OPENMP, CERES_USE_CXX11_THREADS or CERES_NO_THREADS must be defined.
+# error One of CERES_USE_OPENMP, CERES_USE_CXX_THREADS or CERES_NO_THREADS must be defined.
#endif
// CERES_NO_SPARSE should be automatically defined by config.h if Ceres was
// compiled without any sparse back-end. Verify that it has not subsequently
// been inconsistently redefined.
#if defined(CERES_NO_SPARSE)
-# if !defined(CERES_NO_SUITESPARSE)
-# error CERES_NO_SPARSE requires CERES_NO_SUITESPARSE.
-# endif
-# if !defined(CERES_NO_CXSPARSE)
-# error CERES_NO_SPARSE requires CERES_NO_CXSPARSE
-# endif
-# if !defined(CERES_NO_ACCELERATE_SPARSE)
-# error CERES_NO_SPARSE requires CERES_NO_ACCELERATE_SPARSE
-# endif
-# if defined(CERES_USE_EIGEN_SPARSE)
-# error CERES_NO_SPARSE requires !CERES_USE_EIGEN_SPARSE
-# endif
+#if !defined(CERES_NO_SUITESPARSE)
+#error CERES_NO_SPARSE requires CERES_NO_SUITESPARSE.
+#endif
+#if !defined(CERES_NO_CXSPARSE)
+#error CERES_NO_SPARSE requires CERES_NO_CXSPARSE
+#endif
+#if !defined(CERES_NO_ACCELERATE_SPARSE)
+#error CERES_NO_SPARSE requires CERES_NO_ACCELERATE_SPARSE
+#endif
+#if defined(CERES_USE_EIGEN_SPARSE)
+#error CERES_NO_SPARSE requires !CERES_USE_EIGEN_SPARSE
+#endif
#endif
// A macro to signal which functions and classes are exported when
-// building a DLL with MSVC.
-//
-// Note that the ordering here is important, CERES_BUILDING_SHARED_LIBRARY
-// is only defined locally when Ceres is compiled, it is never exported to
-// users. However, in order that we do not have to configure config.h
-// separately for building vs installing, if we are using MSVC and building
-// a shared library, then both CERES_BUILDING_SHARED_LIBRARY and
-// CERES_USING_SHARED_LIBRARY will be defined when Ceres is compiled.
-// Hence it is important that the check for CERES_BUILDING_SHARED_LIBRARY
-// happens first.
-#if defined(_MSC_VER) && defined(CERES_BUILDING_SHARED_LIBRARY)
-# define CERES_EXPORT __declspec(dllexport)
-#elif defined(_MSC_VER) && defined(CERES_USING_SHARED_LIBRARY)
-# define CERES_EXPORT __declspec(dllimport)
+// building a shared library.
+#if defined(_MSC_VER)
+#define CERES_API_SHARED_IMPORT __declspec(dllimport)
+#define CERES_API_SHARED_EXPORT __declspec(dllexport)
+#elif defined(__GNUC__)
+#define CERES_API_SHARED_IMPORT __attribute__((visibility("default")))
+#define CERES_API_SHARED_EXPORT __attribute__((visibility("default")))
#else
-# define CERES_EXPORT
+#define CERES_API_SHARED_IMPORT
+#define CERES_API_SHARED_EXPORT
+#endif
+
+// CERES_BUILDING_SHARED_LIBRARY is only defined locally when Ceres itself is
+// compiled as a shared library, it is never exported to users. In order that
+// we do not have to configure config.h separately when building Ceres as either
+// a static or dynamic library, we define both CERES_USING_SHARED_LIBRARY and
+// CERES_BUILDING_SHARED_LIBRARY when building as a shared library.
+#if defined(CERES_USING_SHARED_LIBRARY)
+#if defined(CERES_BUILDING_SHARED_LIBRARY)
+// Compiling Ceres itself as a shared library.
+#define CERES_EXPORT CERES_API_SHARED_EXPORT
+#else
+// Using Ceres as a shared library.
+#define CERES_EXPORT CERES_API_SHARED_IMPORT
+#endif
+#else
+// Ceres was compiled as a static library, export everything.
+#define CERES_EXPORT
+#endif
+
+// Unit tests reach in and test internal functionality so we need a way to make
+// those symbols visible
+#ifdef CERES_EXPORT_INTERNAL_SYMBOLS
+#define CERES_EXPORT_INTERNAL CERES_EXPORT
+#else
+#define CERES_EXPORT_INTERNAL
#endif
#endif // CERES_PUBLIC_INTERNAL_PORT_H_
diff --git a/include/ceres/internal/reenable_warnings.h b/include/ceres/internal/reenable_warnings.h
index 7e41025..2c5db06 100644
--- a/include/ceres/internal/reenable_warnings.h
+++ b/include/ceres/internal/reenable_warnings.h
@@ -32,7 +32,7 @@
#undef CERES_WARNINGS_DISABLED
#ifdef _MSC_VER
-#pragma warning( pop )
+#pragma warning(pop)
#endif
#endif // CERES_WARNINGS_DISABLED
diff --git a/include/ceres/internal/variadic_evaluate.h b/include/ceres/internal/variadic_evaluate.h
index 26428d0..47ff6b1 100644
--- a/include/ceres/internal/variadic_evaluate.h
+++ b/include/ceres/internal/variadic_evaluate.h
@@ -36,6 +36,7 @@
#include <stddef.h>
#include <type_traits>
+#include <utility>
#include "ceres/cost_function.h"
#include "ceres/internal/parameter_dims.h"
@@ -45,9 +46,11 @@
// For fixed size cost functors
template <typename Functor, typename T, int... Indices>
-inline bool VariadicEvaluateImpl(const Functor& functor, T const* const* input,
- T* output, std::false_type /*is_dynamic*/,
- integer_sequence<int, Indices...>) {
+inline bool VariadicEvaluateImpl(const Functor& functor,
+ T const* const* input,
+ T* output,
+ std::false_type /*is_dynamic*/,
+ std::integer_sequence<int, Indices...>) {
static_assert(sizeof...(Indices),
"Invalid number of parameter blocks. At least one parameter "
"block must be specified.");
@@ -56,26 +59,31 @@
// For dynamic sized cost functors
template <typename Functor, typename T>
-inline bool VariadicEvaluateImpl(const Functor& functor, T const* const* input,
- T* output, std::true_type /*is_dynamic*/,
- integer_sequence<int>) {
+inline bool VariadicEvaluateImpl(const Functor& functor,
+ T const* const* input,
+ T* output,
+ std::true_type /*is_dynamic*/,
+ std::integer_sequence<int>) {
return functor(input, output);
}
// For ceres cost functors (not ceres::CostFunction)
template <typename ParameterDims, typename Functor, typename T>
-inline bool VariadicEvaluateImpl(const Functor& functor, T const* const* input,
- T* output, const void* /* NOT USED */) {
+inline bool VariadicEvaluateImpl(const Functor& functor,
+ T const* const* input,
+ T* output,
+ const void* /* NOT USED */) {
using ParameterBlockIndices =
- make_integer_sequence<int, ParameterDims::kNumParameterBlocks>;
+ std::make_integer_sequence<int, ParameterDims::kNumParameterBlocks>;
using IsDynamic = std::integral_constant<bool, ParameterDims::kIsDynamic>;
- return VariadicEvaluateImpl(functor, input, output, IsDynamic(),
- ParameterBlockIndices());
+ return VariadicEvaluateImpl(
+ functor, input, output, IsDynamic(), ParameterBlockIndices());
}
// For ceres::CostFunction
template <typename ParameterDims, typename Functor, typename T>
-inline bool VariadicEvaluateImpl(const Functor& functor, T const* const* input,
+inline bool VariadicEvaluateImpl(const Functor& functor,
+ T const* const* input,
T* output,
const CostFunction* /* NOT USED */) {
return functor.Evaluate(input, output, nullptr);
@@ -94,7 +102,8 @@
// blocks. The signature of the functor must have the following signature
// 'bool()(const T* i_1, const T* i_2, ... const T* i_n, T* output)'.
template <typename ParameterDims, typename Functor, typename T>
-inline bool VariadicEvaluate(const Functor& functor, T const* const* input,
+inline bool VariadicEvaluate(const Functor& functor,
+ T const* const* input,
T* output) {
return VariadicEvaluateImpl<ParameterDims>(functor, input, output, &functor);
}
diff --git a/include/ceres/iteration_callback.h b/include/ceres/iteration_callback.h
index bd1e782..4507fdf 100644
--- a/include/ceres/iteration_callback.h
+++ b/include/ceres/iteration_callback.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -35,42 +35,22 @@
#ifndef CERES_PUBLIC_ITERATION_CALLBACK_H_
#define CERES_PUBLIC_ITERATION_CALLBACK_H_
-#include "ceres/types.h"
#include "ceres/internal/disable_warnings.h"
+#include "ceres/types.h"
namespace ceres {
// This struct describes the state of the optimizer after each
// iteration of the minimization.
struct CERES_EXPORT IterationSummary {
- IterationSummary()
- : iteration(0),
- step_is_valid(false),
- step_is_nonmonotonic(false),
- step_is_successful(false),
- cost(0.0),
- cost_change(0.0),
- gradient_max_norm(0.0),
- gradient_norm(0.0),
- step_norm(0.0),
- eta(0.0),
- step_size(0.0),
- line_search_function_evaluations(0),
- line_search_gradient_evaluations(0),
- line_search_iterations(0),
- linear_solver_iterations(0),
- iteration_time_in_seconds(0.0),
- step_solver_time_in_seconds(0.0),
- cumulative_time_in_seconds(0.0) {}
-
// Current iteration number.
- int iteration;
+ int iteration = 0;
// Step was numerically valid, i.e., all values are finite and the
// step reduces the value of the linearized model.
//
// Note: step_is_valid is always true when iteration = 0.
- bool step_is_valid;
+ bool step_is_valid = false;
// Step did not reduce the value of the objective function
// sufficiently, but it was accepted because of the relaxed
@@ -78,7 +58,7 @@
// algorithm.
//
// Note: step_is_nonmonotonic is always false when iteration = 0;
- bool step_is_nonmonotonic;
+ bool step_is_nonmonotonic = false;
// Whether or not the minimizer accepted this step or not. If the
// ordinary trust region algorithm is used, this means that the
@@ -90,68 +70,68 @@
// step and the step is declared successful.
//
// Note: step_is_successful is always true when iteration = 0.
- bool step_is_successful;
+ bool step_is_successful = false;
// Value of the objective function.
- double cost;
+ double cost = 0.0;
// Change in the value of the objective function in this
// iteration. This can be positive or negative.
- double cost_change;
+ double cost_change = 0.0;
// Infinity norm of the gradient vector.
- double gradient_max_norm;
+ double gradient_max_norm = 0.0;
// 2-norm of the gradient vector.
- double gradient_norm;
+ double gradient_norm = 0.0;
// 2-norm of the size of the step computed by the optimization
// algorithm.
- double step_norm;
+ double step_norm = 0.0;
// For trust region algorithms, the ratio of the actual change in
// cost and the change in the cost of the linearized approximation.
- double relative_decrease;
+ double relative_decrease = 0.0;
// Size of the trust region at the end of the current iteration. For
// the Levenberg-Marquardt algorithm, the regularization parameter
// mu = 1.0 / trust_region_radius.
- double trust_region_radius;
+ double trust_region_radius = 0.0;
// For the inexact step Levenberg-Marquardt algorithm, this is the
// relative accuracy with which the Newton(LM) step is solved. This
// number affects only the iterative solvers capable of solving
// linear systems inexactly. Factorization-based exact solvers
// ignore it.
- double eta;
+ double eta = 0.0;
// Step sized computed by the line search algorithm.
- double step_size;
+ double step_size = 0.0;
// Number of function value evaluations used by the line search algorithm.
- int line_search_function_evaluations;
+ int line_search_function_evaluations = 0;
// Number of function gradient evaluations used by the line search algorithm.
- int line_search_gradient_evaluations;
+ int line_search_gradient_evaluations = 0;
// Number of iterations taken by the line search algorithm.
- int line_search_iterations;
+ int line_search_iterations = 0;
// Number of iterations taken by the linear solver to solve for the
// Newton step.
- int linear_solver_iterations;
+ int linear_solver_iterations = 0;
// All times reported below are wall times.
// Time (in seconds) spent inside the minimizer loop in the current
// iteration.
- double iteration_time_in_seconds;
+ double iteration_time_in_seconds = 0.0;
// Time (in seconds) spent inside the trust region step solver.
- double step_solver_time_in_seconds;
+ double step_solver_time_in_seconds = 0.0;
// Time (in seconds) since the user called Solve().
- double cumulative_time_in_seconds;
+ double cumulative_time_in_seconds = 0.0;
};
// Interface for specifying callbacks that are executed at the end of
diff --git a/include/ceres/jet.h b/include/ceres/jet.h
index 2b54064..da49f32 100644
--- a/include/ceres/jet.h
+++ b/include/ceres/jet.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -178,20 +178,18 @@
// (where T is a Jet<T, N>). This usually only happens in opt mode. Note that
// the C++ standard mandates that e.g. default constructed doubles are
// initialized to 0.0; see sections 8.5 of the C++03 standard.
- Jet() : a() {
- v.setZero();
- }
+ Jet() : a() { v.setConstant(Scalar()); }
// Constructor from scalar: a + 0.
explicit Jet(const T& value) {
a = value;
- v.setZero();
+ v.setConstant(Scalar());
}
// Constructor from scalar plus variable: a + t_i.
Jet(const T& value, int k) {
a = value;
- v.setZero();
+ v.setConstant(Scalar());
v[k] = T(1.0);
}
@@ -199,28 +197,27 @@
// The use of Eigen::DenseBase allows Eigen expressions
// to be passed in without being fully evaluated until
// they are assigned to v
- template<typename Derived>
- EIGEN_STRONG_INLINE Jet(const T& a, const Eigen::DenseBase<Derived> &v)
- : a(a), v(v) {
- }
+ template <typename Derived>
+ EIGEN_STRONG_INLINE Jet(const T& a, const Eigen::DenseBase<Derived>& v)
+ : a(a), v(v) {}
// Compound operators
- Jet<T, N>& operator+=(const Jet<T, N> &y) {
+ Jet<T, N>& operator+=(const Jet<T, N>& y) {
*this = *this + y;
return *this;
}
- Jet<T, N>& operator-=(const Jet<T, N> &y) {
+ Jet<T, N>& operator-=(const Jet<T, N>& y) {
*this = *this - y;
return *this;
}
- Jet<T, N>& operator*=(const Jet<T, N> &y) {
+ Jet<T, N>& operator*=(const Jet<T, N>& y) {
*this = *this * y;
return *this;
}
- Jet<T, N>& operator/=(const Jet<T, N> &y) {
+ Jet<T, N>& operator/=(const Jet<T, N>& y) {
*this = *this / y;
return *this;
}
@@ -250,66 +247,16 @@
T a;
// The infinitesimal part.
- //
- // We allocate Jets on the stack and other places they might not be aligned
- // to X(=16 [SSE], 32 [AVX] etc)-byte boundaries, which would prevent the safe
- // use of vectorisation. If we have C++11, we can specify the alignment.
- // However, the standard gives wide latitude as to what alignments are valid,
- // and it might be that the maximum supported alignment *guaranteed* to be
- // supported is < 16, in which case we do not specify an alignment, as this
- // implies the host is not a modern x86 machine. If using < C++11, we cannot
- // specify alignment.
+ Eigen::Matrix<T, N, 1> v;
-#if defined(EIGEN_DONT_VECTORIZE)
- Eigen::Matrix<T, N, 1, Eigen::DontAlign> v;
-#else
- // Enable vectorisation iff the maximum supported scalar alignment is >=
- // 16 bytes, as this is the minimum required by Eigen for any vectorisation.
- //
- // NOTE: It might be the case that we could get >= 16-byte alignment even if
- // max_align_t < 16. However we can't guarantee that this
- // would happen (and it should not for any modern x86 machine) and if it
- // didn't, we could get misaligned Jets.
- static constexpr int kAlignOrNot =
- // Work around a GCC 4.8 bug
- // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=56019) where
- // std::max_align_t is misplaced.
-#if defined (__GNUC__) && __GNUC__ == 4 && __GNUC_MINOR__ == 8
- alignof(::max_align_t) >= 16
-#else
- alignof(std::max_align_t) >= 16
-#endif
- ? Eigen::AutoAlign : Eigen::DontAlign;
-
-#if defined(EIGEN_MAX_ALIGN_BYTES)
- // Eigen >= 3.3 supports AVX & FMA instructions that require 32-byte alignment
- // (greater for AVX512). Rather than duplicating the detection logic, use
- // Eigen's macro for the alignment size.
- //
- // NOTE: EIGEN_MAX_ALIGN_BYTES can be > 16 (e.g. 32 for AVX), even though
- // kMaxAlignBytes will max out at 16. We are therefore relying on
- // Eigen's detection logic to ensure that this does not result in
- // misaligned Jets.
-#define CERES_JET_ALIGN_BYTES EIGEN_MAX_ALIGN_BYTES
-#else
- // Eigen < 3.3 only supported 16-byte alignment.
-#define CERES_JET_ALIGN_BYTES 16
-#endif
-
- // Default to the native alignment if 16-byte alignment is not guaranteed to
- // be supported. We cannot use alignof(T) as if we do, GCC 4.8 complains that
- // the alignment 'is not an integer constant', although Clang accepts it.
- static constexpr size_t kAlignment = kAlignOrNot == Eigen::AutoAlign
- ? CERES_JET_ALIGN_BYTES : alignof(double);
-
-#undef CERES_JET_ALIGN_BYTES
- alignas(kAlignment) Eigen::Matrix<T, N, 1, kAlignOrNot> v;
-#endif
+ // This struct needs to have an Eigen aligned operator new as it contains
+ // fixed-size Eigen types.
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW
};
// Unary +
-template<typename T, int N> inline
-Jet<T, N> const& operator+(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> const& operator+(const Jet<T, N>& f) {
return f;
}
@@ -317,72 +264,68 @@
// see if it causes a performance increase.
// Unary -
-template<typename T, int N> inline
-Jet<T, N> operator-(const Jet<T, N>&f) {
+template <typename T, int N>
+inline Jet<T, N> operator-(const Jet<T, N>& f) {
return Jet<T, N>(-f.a, -f.v);
}
// Binary +
-template<typename T, int N> inline
-Jet<T, N> operator+(const Jet<T, N>& f,
- const Jet<T, N>& g) {
+template <typename T, int N>
+inline Jet<T, N> operator+(const Jet<T, N>& f, const Jet<T, N>& g) {
return Jet<T, N>(f.a + g.a, f.v + g.v);
}
// Binary + with a scalar: x + s
-template<typename T, int N> inline
-Jet<T, N> operator+(const Jet<T, N>& f, T s) {
+template <typename T, int N>
+inline Jet<T, N> operator+(const Jet<T, N>& f, T s) {
return Jet<T, N>(f.a + s, f.v);
}
// Binary + with a scalar: s + x
-template<typename T, int N> inline
-Jet<T, N> operator+(T s, const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> operator+(T s, const Jet<T, N>& f) {
return Jet<T, N>(f.a + s, f.v);
}
// Binary -
-template<typename T, int N> inline
-Jet<T, N> operator-(const Jet<T, N>& f,
- const Jet<T, N>& g) {
+template <typename T, int N>
+inline Jet<T, N> operator-(const Jet<T, N>& f, const Jet<T, N>& g) {
return Jet<T, N>(f.a - g.a, f.v - g.v);
}
// Binary - with a scalar: x - s
-template<typename T, int N> inline
-Jet<T, N> operator-(const Jet<T, N>& f, T s) {
+template <typename T, int N>
+inline Jet<T, N> operator-(const Jet<T, N>& f, T s) {
return Jet<T, N>(f.a - s, f.v);
}
// Binary - with a scalar: s - x
-template<typename T, int N> inline
-Jet<T, N> operator-(T s, const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> operator-(T s, const Jet<T, N>& f) {
return Jet<T, N>(s - f.a, -f.v);
}
// Binary *
-template<typename T, int N> inline
-Jet<T, N> operator*(const Jet<T, N>& f,
- const Jet<T, N>& g) {
+template <typename T, int N>
+inline Jet<T, N> operator*(const Jet<T, N>& f, const Jet<T, N>& g) {
return Jet<T, N>(f.a * g.a, f.a * g.v + f.v * g.a);
}
// Binary * with a scalar: x * s
-template<typename T, int N> inline
-Jet<T, N> operator*(const Jet<T, N>& f, T s) {
+template <typename T, int N>
+inline Jet<T, N> operator*(const Jet<T, N>& f, T s) {
return Jet<T, N>(f.a * s, f.v * s);
}
// Binary * with a scalar: s * x
-template<typename T, int N> inline
-Jet<T, N> operator*(T s, const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> operator*(T s, const Jet<T, N>& f) {
return Jet<T, N>(f.a * s, f.v * s);
}
// Binary /
-template<typename T, int N> inline
-Jet<T, N> operator/(const Jet<T, N>& f,
- const Jet<T, N>& g) {
+template <typename T, int N>
+inline Jet<T, N> operator/(const Jet<T, N>& f, const Jet<T, N>& g) {
// This uses:
//
// a + u (a + u)(b - v) (a + u)(b - v)
@@ -396,39 +339,39 @@
}
// Binary / with a scalar: s / x
-template<typename T, int N> inline
-Jet<T, N> operator/(T s, const Jet<T, N>& g) {
+template <typename T, int N>
+inline Jet<T, N> operator/(T s, const Jet<T, N>& g) {
const T minus_s_g_a_inverse2 = -s / (g.a * g.a);
return Jet<T, N>(s / g.a, g.v * minus_s_g_a_inverse2);
}
// Binary / with a scalar: x / s
-template<typename T, int N> inline
-Jet<T, N> operator/(const Jet<T, N>& f, T s) {
+template <typename T, int N>
+inline Jet<T, N> operator/(const Jet<T, N>& f, T s) {
const T s_inverse = T(1.0) / s;
return Jet<T, N>(f.a * s_inverse, f.v * s_inverse);
}
// Binary comparison operators for both scalars and jets.
-#define CERES_DEFINE_JET_COMPARISON_OPERATOR(op) \
-template<typename T, int N> inline \
-bool operator op(const Jet<T, N>& f, const Jet<T, N>& g) { \
- return f.a op g.a; \
-} \
-template<typename T, int N> inline \
-bool operator op(const T& s, const Jet<T, N>& g) { \
- return s op g.a; \
-} \
-template<typename T, int N> inline \
-bool operator op(const Jet<T, N>& f, const T& s) { \
- return f.a op s; \
-}
-CERES_DEFINE_JET_COMPARISON_OPERATOR( < ) // NOLINT
-CERES_DEFINE_JET_COMPARISON_OPERATOR( <= ) // NOLINT
-CERES_DEFINE_JET_COMPARISON_OPERATOR( > ) // NOLINT
-CERES_DEFINE_JET_COMPARISON_OPERATOR( >= ) // NOLINT
-CERES_DEFINE_JET_COMPARISON_OPERATOR( == ) // NOLINT
-CERES_DEFINE_JET_COMPARISON_OPERATOR( != ) // NOLINT
+#define CERES_DEFINE_JET_COMPARISON_OPERATOR(op) \
+ template <typename T, int N> \
+ inline bool operator op(const Jet<T, N>& f, const Jet<T, N>& g) { \
+ return f.a op g.a; \
+ } \
+ template <typename T, int N> \
+ inline bool operator op(const T& s, const Jet<T, N>& g) { \
+ return s op g.a; \
+ } \
+ template <typename T, int N> \
+ inline bool operator op(const Jet<T, N>& f, const T& s) { \
+ return f.a op s; \
+ }
+CERES_DEFINE_JET_COMPARISON_OPERATOR(<) // NOLINT
+CERES_DEFINE_JET_COMPARISON_OPERATOR(<=) // NOLINT
+CERES_DEFINE_JET_COMPARISON_OPERATOR(>) // NOLINT
+CERES_DEFINE_JET_COMPARISON_OPERATOR(>=) // NOLINT
+CERES_DEFINE_JET_COMPARISON_OPERATOR(==) // NOLINT
+CERES_DEFINE_JET_COMPARISON_OPERATOR(!=) // NOLINT
#undef CERES_DEFINE_JET_COMPARISON_OPERATOR
// Pull some functions from namespace std.
@@ -445,6 +388,8 @@
using std::ceil;
using std::cos;
using std::cosh;
+using std::erf;
+using std::erfc;
using std::exp;
using std::exp2;
using std::floor;
@@ -465,97 +410,99 @@
using std::tanh;
// Legacy names from pre-C++11 days.
-inline bool IsFinite (double x) { return std::isfinite(x); }
+// clang-format off
+inline bool IsFinite(double x) { return std::isfinite(x); }
inline bool IsInfinite(double x) { return std::isinf(x); }
-inline bool IsNaN (double x) { return std::isnan(x); }
-inline bool IsNormal (double x) { return std::isnormal(x); }
+inline bool IsNaN(double x) { return std::isnan(x); }
+inline bool IsNormal(double x) { return std::isnormal(x); }
+// clang-format on
// In general, f(a + h) ~= f(a) + f'(a) h, via the chain rule.
// abs(x + h) ~= x + h or -(x + h)
-template <typename T, int N> inline
-Jet<T, N> abs(const Jet<T, N>& f) {
- return f.a < T(0.0) ? -f : f;
+template <typename T, int N>
+inline Jet<T, N> abs(const Jet<T, N>& f) {
+ return (f.a < T(0.0) ? -f : f);
}
// log(a + h) ~= log(a) + h / a
-template <typename T, int N> inline
-Jet<T, N> log(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> log(const Jet<T, N>& f) {
const T a_inverse = T(1.0) / f.a;
return Jet<T, N>(log(f.a), f.v * a_inverse);
}
// exp(a + h) ~= exp(a) + exp(a) h
-template <typename T, int N> inline
-Jet<T, N> exp(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> exp(const Jet<T, N>& f) {
const T tmp = exp(f.a);
return Jet<T, N>(tmp, tmp * f.v);
}
// sqrt(a + h) ~= sqrt(a) + h / (2 sqrt(a))
-template <typename T, int N> inline
-Jet<T, N> sqrt(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> sqrt(const Jet<T, N>& f) {
const T tmp = sqrt(f.a);
const T two_a_inverse = T(1.0) / (T(2.0) * tmp);
return Jet<T, N>(tmp, f.v * two_a_inverse);
}
// cos(a + h) ~= cos(a) - sin(a) h
-template <typename T, int N> inline
-Jet<T, N> cos(const Jet<T, N>& f) {
- return Jet<T, N>(cos(f.a), - sin(f.a) * f.v);
+template <typename T, int N>
+inline Jet<T, N> cos(const Jet<T, N>& f) {
+ return Jet<T, N>(cos(f.a), -sin(f.a) * f.v);
}
// acos(a + h) ~= acos(a) - 1 / sqrt(1 - a^2) h
-template <typename T, int N> inline
-Jet<T, N> acos(const Jet<T, N>& f) {
- const T tmp = - T(1.0) / sqrt(T(1.0) - f.a * f.a);
+template <typename T, int N>
+inline Jet<T, N> acos(const Jet<T, N>& f) {
+ const T tmp = -T(1.0) / sqrt(T(1.0) - f.a * f.a);
return Jet<T, N>(acos(f.a), tmp * f.v);
}
// sin(a + h) ~= sin(a) + cos(a) h
-template <typename T, int N> inline
-Jet<T, N> sin(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> sin(const Jet<T, N>& f) {
return Jet<T, N>(sin(f.a), cos(f.a) * f.v);
}
// asin(a + h) ~= asin(a) + 1 / sqrt(1 - a^2) h
-template <typename T, int N> inline
-Jet<T, N> asin(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> asin(const Jet<T, N>& f) {
const T tmp = T(1.0) / sqrt(T(1.0) - f.a * f.a);
return Jet<T, N>(asin(f.a), tmp * f.v);
}
// tan(a + h) ~= tan(a) + (1 + tan(a)^2) h
-template <typename T, int N> inline
-Jet<T, N> tan(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> tan(const Jet<T, N>& f) {
const T tan_a = tan(f.a);
const T tmp = T(1.0) + tan_a * tan_a;
return Jet<T, N>(tan_a, tmp * f.v);
}
// atan(a + h) ~= atan(a) + 1 / (1 + a^2) h
-template <typename T, int N> inline
-Jet<T, N> atan(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> atan(const Jet<T, N>& f) {
const T tmp = T(1.0) / (T(1.0) + f.a * f.a);
return Jet<T, N>(atan(f.a), tmp * f.v);
}
// sinh(a + h) ~= sinh(a) + cosh(a) h
-template <typename T, int N> inline
-Jet<T, N> sinh(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> sinh(const Jet<T, N>& f) {
return Jet<T, N>(sinh(f.a), cosh(f.a) * f.v);
}
// cosh(a + h) ~= cosh(a) + sinh(a) h
-template <typename T, int N> inline
-Jet<T, N> cosh(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> cosh(const Jet<T, N>& f) {
return Jet<T, N>(cosh(f.a), sinh(f.a) * f.v);
}
// tanh(a + h) ~= tanh(a) + (1 - tanh(a)^2) h
-template <typename T, int N> inline
-Jet<T, N> tanh(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> tanh(const Jet<T, N>& f) {
const T tanh_a = tanh(f.a);
const T tmp = T(1.0) - tanh_a * tanh_a;
return Jet<T, N>(tanh_a, tmp * f.v);
@@ -565,8 +512,8 @@
// result in a zero derivative which provides no information to the solver.
//
// floor(a + h) ~= floor(a) + 0
-template <typename T, int N> inline
-Jet<T, N> floor(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> floor(const Jet<T, N>& f) {
return Jet<T, N>(floor(f.a));
}
@@ -574,31 +521,31 @@
// result in a zero derivative which provides no information to the solver.
//
// ceil(a + h) ~= ceil(a) + 0
-template <typename T, int N> inline
-Jet<T, N> ceil(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> ceil(const Jet<T, N>& f) {
return Jet<T, N>(ceil(f.a));
}
// Some new additions to C++11:
// cbrt(a + h) ~= cbrt(a) + h / (3 a ^ (2/3))
-template <typename T, int N> inline
-Jet<T, N> cbrt(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> cbrt(const Jet<T, N>& f) {
const T derivative = T(1.0) / (T(3.0) * cbrt(f.a * f.a));
return Jet<T, N>(cbrt(f.a), f.v * derivative);
}
// exp2(x + h) = 2^(x+h) ~= 2^x + h*2^x*log(2)
-template <typename T, int N> inline
-Jet<T, N> exp2(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> exp2(const Jet<T, N>& f) {
const T tmp = exp2(f.a);
const T derivative = tmp * log(T(2));
return Jet<T, N>(tmp, f.v * derivative);
}
// log2(x + h) ~= log2(x) + h / (x * log(2))
-template <typename T, int N> inline
-Jet<T, N> log2(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> log2(const Jet<T, N>& f) {
const T derivative = T(1.0) / (f.a * log(T(2)));
return Jet<T, N>(log2(f.a), f.v * derivative);
}
@@ -607,8 +554,8 @@
// but acts to prevent underflow/overflow for small/large x/y.
// Note that the function is non-smooth at x=y=0,
// so the derivative is undefined there.
-template <typename T, int N> inline
-Jet<T, N> hypot(const Jet<T, N>& x, const Jet<T, N>& y) {
+template <typename T, int N>
+inline Jet<T, N> hypot(const Jet<T, N>& x, const Jet<T, N>& y) {
// d/da sqrt(a) = 0.5 / sqrt(a)
// d/dx x^2 + y^2 = 2x
// So by the chain rule:
@@ -618,16 +565,31 @@
return Jet<T, N>(tmp, x.a / tmp * x.v + y.a / tmp * y.v);
}
-template <typename T, int N> inline
-const Jet<T, N>& fmax(const Jet<T, N>& x, const Jet<T, N>& y) {
+template <typename T, int N>
+inline Jet<T, N> fmax(const Jet<T, N>& x, const Jet<T, N>& y) {
return x < y ? y : x;
}
-template <typename T, int N> inline
-const Jet<T, N>& fmin(const Jet<T, N>& x, const Jet<T, N>& y) {
+template <typename T, int N>
+inline Jet<T, N> fmin(const Jet<T, N>& x, const Jet<T, N>& y) {
return y < x ? y : x;
}
+// erf is defined as an integral that cannot be expressed analyticaly
+// however, the derivative is trivial to compute
+// erf(x + h) = erf(x) + h * 2*exp(-x^2)/sqrt(pi)
+template <typename T, int N>
+inline Jet<T, N> erf(const Jet<T, N>& x) {
+ return Jet<T, N>(erf(x.a), x.v * M_2_SQRTPI * exp(-x.a * x.a));
+}
+
+// erfc(x) = 1-erf(x)
+// erfc(x + h) = erfc(x) + h * (-2*exp(-x^2)/sqrt(pi))
+template <typename T, int N>
+inline Jet<T, N> erfc(const Jet<T, N>& x) {
+ return Jet<T, N>(erfc(x.a), -x.v * M_2_SQRTPI * exp(-x.a * x.a));
+}
+
// Bessel functions of the first kind with integer order equal to 0, 1, n.
//
// Microsoft has deprecated the j[0,1,n]() POSIX Bessel functions in favour of
@@ -664,32 +626,32 @@
// See formula http://dlmf.nist.gov/10.6#E3
// j0(a + h) ~= j0(a) - j1(a) h
-template <typename T, int N> inline
-Jet<T, N> BesselJ0(const Jet<T, N>& f) {
- return Jet<T, N>(BesselJ0(f.a),
- -BesselJ1(f.a) * f.v);
+template <typename T, int N>
+inline Jet<T, N> BesselJ0(const Jet<T, N>& f) {
+ return Jet<T, N>(BesselJ0(f.a), -BesselJ1(f.a) * f.v);
}
// See formula http://dlmf.nist.gov/10.6#E1
// j1(a + h) ~= j1(a) + 0.5 ( j0(a) - j2(a) ) h
-template <typename T, int N> inline
-Jet<T, N> BesselJ1(const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> BesselJ1(const Jet<T, N>& f) {
return Jet<T, N>(BesselJ1(f.a),
T(0.5) * (BesselJ0(f.a) - BesselJn(2, f.a)) * f.v);
}
// See formula http://dlmf.nist.gov/10.6#E1
// j_n(a + h) ~= j_n(a) + 0.5 ( j_{n-1}(a) - j_{n+1}(a) ) h
-template <typename T, int N> inline
-Jet<T, N> BesselJn(int n, const Jet<T, N>& f) {
- return Jet<T, N>(BesselJn(n, f.a),
- T(0.5) * (BesselJn(n - 1, f.a) - BesselJn(n + 1, f.a)) * f.v);
+template <typename T, int N>
+inline Jet<T, N> BesselJn(int n, const Jet<T, N>& f) {
+ return Jet<T, N>(
+ BesselJn(n, f.a),
+ T(0.5) * (BesselJn(n - 1, f.a) - BesselJn(n + 1, f.a)) * f.v);
}
// Jet Classification. It is not clear what the appropriate semantics are for
-// these classifications. This picks that std::isfinite and std::isnormal are "all"
-// operations, i.e. all elements of the jet must be finite for the jet itself
-// to be finite (or normal). For IsNaN and IsInfinite, the answer is less
+// these classifications. This picks that std::isfinite and std::isnormal are
+// "all" operations, i.e. all elements of the jet must be finite for the jet
+// itself to be finite (or normal). For IsNaN and IsInfinite, the answer is less
// clear. This takes a "any" approach for IsNaN and IsInfinite such that if any
// part of a jet is nan or inf, then the entire jet is nan or inf. This leads
// to strange situations like a jet can be both IsInfinite and IsNaN, but in
@@ -697,60 +659,45 @@
// derivatives are sane.
// The jet is finite if all parts of the jet are finite.
-template <typename T, int N> inline
-bool isfinite(const Jet<T, N>& f) {
- if (!std::isfinite(f.a)) {
- return false;
- }
+template <typename T, int N>
+inline bool isfinite(const Jet<T, N>& f) {
+ // Branchless implementation. This is more efficient for the false-case and
+ // works with the codegen system.
+ auto result = isfinite(f.a);
for (int i = 0; i < N; ++i) {
- if (!std::isfinite(f.v[i])) {
- return false;
- }
+ result = result & isfinite(f.v[i]);
}
- return true;
+ return result;
}
// The jet is infinite if any part of the Jet is infinite.
-template <typename T, int N> inline
-bool isinf(const Jet<T, N>& f) {
- if (std::isinf(f.a)) {
- return true;
- }
+template <typename T, int N>
+inline bool isinf(const Jet<T, N>& f) {
+ auto result = isinf(f.a);
for (int i = 0; i < N; ++i) {
- if (std::isinf(f.v[i])) {
- return true;
- }
+ result = result | isinf(f.v[i]);
}
- return false;
+ return result;
}
-
// The jet is NaN if any part of the jet is NaN.
-template <typename T, int N> inline
-bool isnan(const Jet<T, N>& f) {
- if (std::isnan(f.a)) {
- return true;
- }
+template <typename T, int N>
+inline bool isnan(const Jet<T, N>& f) {
+ auto result = isnan(f.a);
for (int i = 0; i < N; ++i) {
- if (std::isnan(f.v[i])) {
- return true;
- }
+ result = result | isnan(f.v[i]);
}
- return false;
+ return result;
}
// The jet is normal if all parts of the jet are normal.
-template <typename T, int N> inline
-bool isnormal(const Jet<T, N>& f) {
- if (!std::isnormal(f.a)) {
- return false;
- }
+template <typename T, int N>
+inline bool isnormal(const Jet<T, N>& f) {
+ auto result = isnormal(f.a);
for (int i = 0; i < N; ++i) {
- if (!std::isnormal(f.v[i])) {
- return false;
- }
+ result = result & isnormal(f.v[i]);
}
- return true;
+ return result;
}
// Legacy functions from the pre-C++11 days.
@@ -770,8 +717,8 @@
}
// The jet is infinite if any part of the jet is infinite.
-template <typename T, int N> inline
-bool IsInfinite(const Jet<T, N>& f) {
+template <typename T, int N>
+inline bool IsInfinite(const Jet<T, N>& f) {
return isinf(f);
}
@@ -779,22 +726,21 @@
//
// In words: the rate of change of theta is 1/r times the rate of
// change of (x, y) in the positive angular direction.
-template <typename T, int N> inline
-Jet<T, N> atan2(const Jet<T, N>& g, const Jet<T, N>& f) {
+template <typename T, int N>
+inline Jet<T, N> atan2(const Jet<T, N>& g, const Jet<T, N>& f) {
// Note order of arguments:
//
// f = a + da
// g = b + db
T const tmp = T(1.0) / (f.a * f.a + g.a * g.a);
- return Jet<T, N>(atan2(g.a, f.a), tmp * (- g.a * f.v + f.a * g.v));
+ return Jet<T, N>(atan2(g.a, f.a), tmp * (-g.a * f.v + f.a * g.v));
}
-
// pow -- base is a differentiable function, exponent is a constant.
// (a+da)^p ~= a^p + p*a^(p-1) da
-template <typename T, int N> inline
-Jet<T, N> pow(const Jet<T, N>& f, double g) {
+template <typename T, int N>
+inline Jet<T, N> pow(const Jet<T, N>& f, double g) {
T const tmp = g * pow(f.a, g - T(1.0));
return Jet<T, N>(pow(f.a, g), tmp * f.v);
}
@@ -810,26 +756,30 @@
// 3. For f < 0 and integer g we have: (f)^(g + dg) ~= f^g but if dg
// != 0, the derivatives are not defined and we return NaN.
-template <typename T, int N> inline
-Jet<T, N> pow(double f, const Jet<T, N>& g) {
- if (f == 0 && g.a > 0) {
+template <typename T, int N>
+inline Jet<T, N> pow(T f, const Jet<T, N>& g) {
+ Jet<T, N> result;
+
+ if (f == T(0) && g.a > T(0)) {
// Handle case 2.
- return Jet<T, N>(T(0.0));
- }
- if (f < 0 && g.a == floor(g.a)) {
- // Handle case 3.
- Jet<T, N> ret(pow(f, g.a));
- for (int i = 0; i < N; i++) {
- if (g.v[i] != T(0.0)) {
- // Return a NaN when g.v != 0.
- ret.v[i] = std::numeric_limits<T>::quiet_NaN();
+ result = Jet<T, N>(T(0.0));
+ } else {
+ if (f < 0 && g.a == floor(g.a)) { // Handle case 3.
+ result = Jet<T, N>(pow(f, g.a));
+ for (int i = 0; i < N; i++) {
+ if (g.v[i] != T(0.0)) {
+ // Return a NaN when g.v != 0.
+ result.v[i] = std::numeric_limits<T>::quiet_NaN();
+ }
}
+ } else {
+ // Handle case 1.
+ T const tmp = pow(f, g.a);
+ result = Jet<T, N>(tmp, log(f) * tmp * g.v);
}
- return ret;
}
- // Handle case 1.
- T const tmp = pow(f, g.a);
- return Jet<T, N>(tmp, log(f) * tmp * g.v);
+
+ return result;
}
// pow -- both base and exponent are differentiable functions. This has a
@@ -868,41 +818,48 @@
//
// 9. For f < 0, g noninteger: The value and derivatives of f^g are not finite.
-template <typename T, int N> inline
-Jet<T, N> pow(const Jet<T, N>& f, const Jet<T, N>& g) {
- if (f.a == 0 && g.a >= 1) {
+template <typename T, int N>
+inline Jet<T, N> pow(const Jet<T, N>& f, const Jet<T, N>& g) {
+ Jet<T, N> result;
+
+ if (f.a == T(0) && g.a >= T(1)) {
// Handle cases 2 and 3.
- if (g.a > 1) {
- return Jet<T, N>(T(0.0));
+ if (g.a > T(1)) {
+ result = Jet<T, N>(T(0.0));
+ } else {
+ result = f;
}
- return f;
- }
- if (f.a < 0 && g.a == floor(g.a)) {
- // Handle cases 7 and 8.
- T const tmp = g.a * pow(f.a, g.a - T(1.0));
- Jet<T, N> ret(pow(f.a, g.a), tmp * f.v);
- for (int i = 0; i < N; i++) {
- if (g.v[i] != T(0.0)) {
- // Return a NaN when g.v != 0.
- ret.v[i] = std::numeric_limits<T>::quiet_NaN();
+
+ } else {
+ if (f.a < T(0) && g.a == floor(g.a)) {
+ // Handle cases 7 and 8.
+ T const tmp = g.a * pow(f.a, g.a - T(1.0));
+ result = Jet<T, N>(pow(f.a, g.a), tmp * f.v);
+ for (int i = 0; i < N; i++) {
+ if (g.v[i] != T(0.0)) {
+ // Return a NaN when g.v != 0.
+ result.v[i] = T(std::numeric_limits<double>::quiet_NaN());
+ }
}
+ } else {
+ // Handle the remaining cases. For cases 4,5,6,9 we allow the log()
+ // function to generate -HUGE_VAL or NaN, since those cases result in a
+ // nonfinite derivative.
+ T const tmp1 = pow(f.a, g.a);
+ T const tmp2 = g.a * pow(f.a, g.a - T(1.0));
+ T const tmp3 = tmp1 * log(f.a);
+ result = Jet<T, N>(tmp1, tmp2 * f.v + tmp3 * g.v);
}
- return ret;
}
- // Handle the remaining cases. For cases 4,5,6,9 we allow the log() function
- // to generate -HUGE_VAL or NaN, since those cases result in a nonfinite
- // derivative.
- T const tmp1 = pow(f.a, g.a);
- T const tmp2 = g.a * pow(f.a, g.a - T(1.0));
- T const tmp3 = tmp1 * log(f.a);
- return Jet<T, N>(tmp1, tmp2 * f.v + tmp3 * g.v);
+
+ return result;
}
// Note: This has to be in the ceres namespace for argument dependent lookup to
// function correctly. Otherwise statements like CHECK_LE(x, 2.0) fail with
// strange compile errors.
template <typename T, int N>
-inline std::ostream &operator<<(std::ostream &s, const Jet<T, N>& z) {
+inline std::ostream& operator<<(std::ostream& s, const Jet<T, N>& z) {
s << "[" << z.a << " ; ";
for (int i = 0; i < N; ++i) {
s << z.v[i];
@@ -913,14 +870,77 @@
s << "]";
return s;
}
-
} // namespace ceres
+namespace std {
+template <typename T, int N>
+struct numeric_limits<ceres::Jet<T, N>> {
+ static constexpr bool is_specialized = true;
+ static constexpr bool is_signed = std::numeric_limits<T>::is_signed;
+ static constexpr bool is_integer = std::numeric_limits<T>::is_integer;
+ static constexpr bool is_exact = std::numeric_limits<T>::is_exact;
+ static constexpr bool has_infinity = std::numeric_limits<T>::has_infinity;
+ static constexpr bool has_quiet_NaN = std::numeric_limits<T>::has_quiet_NaN;
+ static constexpr bool has_signaling_NaN =
+ std::numeric_limits<T>::has_signaling_NaN;
+ static constexpr bool is_iec559 = std::numeric_limits<T>::is_iec559;
+ static constexpr bool is_bounded = std::numeric_limits<T>::is_bounded;
+ static constexpr bool is_modulo = std::numeric_limits<T>::is_modulo;
+
+ static constexpr std::float_denorm_style has_denorm =
+ std::numeric_limits<T>::has_denorm;
+ static constexpr std::float_round_style round_style =
+ std::numeric_limits<T>::round_style;
+
+ static constexpr int digits = std::numeric_limits<T>::digits;
+ static constexpr int digits10 = std::numeric_limits<T>::digits10;
+ static constexpr int max_digits10 = std::numeric_limits<T>::max_digits10;
+ static constexpr int radix = std::numeric_limits<T>::radix;
+ static constexpr int min_exponent = std::numeric_limits<T>::min_exponent;
+ static constexpr int min_exponent10 = std::numeric_limits<T>::max_exponent10;
+ static constexpr int max_exponent = std::numeric_limits<T>::max_exponent;
+ static constexpr int max_exponent10 = std::numeric_limits<T>::max_exponent10;
+ static constexpr bool traps = std::numeric_limits<T>::traps;
+ static constexpr bool tinyness_before =
+ std::numeric_limits<T>::tinyness_before;
+
+ static constexpr ceres::Jet<T, N> min() noexcept {
+ return ceres::Jet<T, N>(std::numeric_limits<T>::min());
+ }
+ static constexpr ceres::Jet<T, N> lowest() noexcept {
+ return ceres::Jet<T, N>(std::numeric_limits<T>::lowest());
+ }
+ static constexpr ceres::Jet<T, N> epsilon() noexcept {
+ return ceres::Jet<T, N>(std::numeric_limits<T>::epsilon());
+ }
+ static constexpr ceres::Jet<T, N> round_error() noexcept {
+ return ceres::Jet<T, N>(std::numeric_limits<T>::round_error());
+ }
+ static constexpr ceres::Jet<T, N> infinity() noexcept {
+ return ceres::Jet<T, N>(std::numeric_limits<T>::infinity());
+ }
+ static constexpr ceres::Jet<T, N> quiet_NaN() noexcept {
+ return ceres::Jet<T, N>(std::numeric_limits<T>::quiet_NaN());
+ }
+ static constexpr ceres::Jet<T, N> signaling_NaN() noexcept {
+ return ceres::Jet<T, N>(std::numeric_limits<T>::signaling_NaN());
+ }
+ static constexpr ceres::Jet<T, N> denorm_min() noexcept {
+ return ceres::Jet<T, N>(std::numeric_limits<T>::denorm_min());
+ }
+
+ static constexpr ceres::Jet<T, N> max() noexcept {
+ return ceres::Jet<T, N>(std::numeric_limits<T>::max());
+ }
+};
+
+} // namespace std
+
namespace Eigen {
// Creating a specialization of NumTraits enables placing Jet objects inside
// Eigen arrays, getting all the goodness of Eigen combined with autodiff.
-template<typename T, int N>
+template <typename T, int N>
struct NumTraits<ceres::Jet<T, N>> {
typedef ceres::Jet<T, N> Real;
typedef ceres::Jet<T, N> NonInteger;
@@ -949,7 +969,7 @@
RequireInitialization = 1
};
- template<bool Vectorized>
+ template <bool Vectorized>
struct Div {
enum {
#if defined(EIGEN_VECTORIZE_AVX)
@@ -968,7 +988,6 @@
static inline Real lowest() { return Real(-std::numeric_limits<T>::max()); }
};
-#if EIGEN_VERSION_AT_LEAST(3, 3, 0)
// Specifying the return type of binary operations between Jets and scalar types
// allows you to perform matrix/array operations with Eigen matrices and arrays
// such as addition, subtraction, multiplication, and division where one Eigen
@@ -983,7 +1002,6 @@
struct ScalarBinaryOpTraits<T, ceres::Jet<T, N>, BinaryOp> {
typedef ceres::Jet<T, N> ReturnType;
};
-#endif // EIGEN_VERSION_AT_LEAST(3, 3, 0)
} // namespace Eigen
diff --git a/include/ceres/local_parameterization.h b/include/ceres/local_parameterization.h
index 5eed035..ba7579d 100644
--- a/include/ceres/local_parameterization.h
+++ b/include/ceres/local_parameterization.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -32,9 +32,12 @@
#ifndef CERES_PUBLIC_LOCAL_PARAMETERIZATION_H_
#define CERES_PUBLIC_LOCAL_PARAMETERIZATION_H_
+#include <array>
+#include <memory>
#include <vector>
-#include "ceres/internal/port.h"
+
#include "ceres/internal/disable_warnings.h"
+#include "ceres/internal/port.h"
namespace ceres {
@@ -87,8 +90,8 @@
//
// An example that occurs commonly in Structure from Motion problems
// is when camera rotations are parameterized using Quaternion. There,
-// it is useful only make updates orthogonal to that 4-vector defining
-// the quaternion. One way to do this is to let delta be a 3
+// it is useful to only make updates orthogonal to that 4-vector
+// defining the quaternion. One way to do this is to let delta be a 3
// dimensional vector and define Plus to be
//
// Plus(x, delta) = [cos(|delta|), sin(|delta|) delta / |delta|] * x
@@ -96,7 +99,7 @@
// The multiplication between the two 4-vectors on the RHS is the
// standard quaternion product.
//
-// Given g and a point x, optimizing f can now be restated as
+// Given f and a point x, optimizing f can now be restated as
//
// min f(Plus(x, delta))
// delta
@@ -153,17 +156,16 @@
public:
explicit IdentityParameterization(int size);
virtual ~IdentityParameterization() {}
- virtual bool Plus(const double* x,
- const double* delta,
- double* x_plus_delta) const;
- virtual bool ComputeJacobian(const double* x,
- double* jacobian) const;
- virtual bool MultiplyByJacobian(const double* x,
- const int num_cols,
- const double* global_matrix,
- double* local_matrix) const;
- virtual int GlobalSize() const { return size_; }
- virtual int LocalSize() const { return size_; }
+ bool Plus(const double* x,
+ const double* delta,
+ double* x_plus_delta) const override;
+ bool ComputeJacobian(const double* x, double* jacobian) const override;
+ bool MultiplyByJacobian(const double* x,
+ const int num_cols,
+ const double* global_matrix,
+ double* local_matrix) const override;
+ int GlobalSize() const override { return size_; }
+ int LocalSize() const override { return size_; }
private:
const int size_;
@@ -175,19 +177,18 @@
explicit SubsetParameterization(int size,
const std::vector<int>& constant_parameters);
virtual ~SubsetParameterization() {}
- virtual bool Plus(const double* x,
- const double* delta,
- double* x_plus_delta) const;
- virtual bool ComputeJacobian(const double* x,
- double* jacobian) const;
- virtual bool MultiplyByJacobian(const double* x,
- const int num_cols,
- const double* global_matrix,
- double* local_matrix) const;
- virtual int GlobalSize() const {
+ bool Plus(const double* x,
+ const double* delta,
+ double* x_plus_delta) const override;
+ bool ComputeJacobian(const double* x, double* jacobian) const override;
+ bool MultiplyByJacobian(const double* x,
+ const int num_cols,
+ const double* global_matrix,
+ double* local_matrix) const override;
+ int GlobalSize() const override {
return static_cast<int>(constancy_mask_.size());
}
- virtual int LocalSize() const { return local_size_; }
+ int LocalSize() const override { return local_size_; }
private:
const int local_size_;
@@ -201,13 +202,12 @@
class CERES_EXPORT QuaternionParameterization : public LocalParameterization {
public:
virtual ~QuaternionParameterization() {}
- virtual bool Plus(const double* x,
- const double* delta,
- double* x_plus_delta) const;
- virtual bool ComputeJacobian(const double* x,
- double* jacobian) const;
- virtual int GlobalSize() const { return 4; }
- virtual int LocalSize() const { return 3; }
+ bool Plus(const double* x,
+ const double* delta,
+ double* x_plus_delta) const override;
+ bool ComputeJacobian(const double* x, double* jacobian) const override;
+ int GlobalSize() const override { return 4; }
+ int LocalSize() const override { return 3; }
};
// Implements the quaternion local parameterization for Eigen's representation
@@ -225,12 +225,12 @@
: public ceres::LocalParameterization {
public:
virtual ~EigenQuaternionParameterization() {}
- virtual bool Plus(const double* x,
- const double* delta,
- double* x_plus_delta) const;
- virtual bool ComputeJacobian(const double* x, double* jacobian) const;
- virtual int GlobalSize() const { return 4; }
- virtual int LocalSize() const { return 3; }
+ bool Plus(const double* x,
+ const double* delta,
+ double* x_plus_delta) const override;
+ bool ComputeJacobian(const double* x, double* jacobian) const override;
+ int GlobalSize() const override { return 4; }
+ int LocalSize() const override { return 3; }
};
// This provides a parameterization for homogeneous vectors which are commonly
@@ -246,32 +246,55 @@
// remain on the sphere. We assume that the last element of x is the scalar
// component. The size of the homogeneous vector is required to be greater than
// 1.
-class CERES_EXPORT HomogeneousVectorParameterization :
- public LocalParameterization {
+class CERES_EXPORT HomogeneousVectorParameterization
+ : public LocalParameterization {
public:
explicit HomogeneousVectorParameterization(int size);
virtual ~HomogeneousVectorParameterization() {}
- virtual bool Plus(const double* x,
- const double* delta,
- double* x_plus_delta) const;
- virtual bool ComputeJacobian(const double* x,
- double* jacobian) const;
- virtual int GlobalSize() const { return size_; }
- virtual int LocalSize() const { return size_ - 1; }
+ bool Plus(const double* x,
+ const double* delta,
+ double* x_plus_delta) const override;
+ bool ComputeJacobian(const double* x, double* jacobian) const override;
+ int GlobalSize() const override { return size_; }
+ int LocalSize() const override { return size_ - 1; }
private:
const int size_;
};
+// This provides a parameterization for lines, where the line is
+// over-parameterized by an origin point and a direction vector. So the
+// parameter vector size needs to be two times the ambient space dimension,
+// where the first half is interpreted as the origin point and the second half
+// as the direction.
+//
+// The plus operator for the line direction is the same as for the
+// HomogeneousVectorParameterization. The update of the origin point is
+// perpendicular to the line direction before the update.
+//
+// This local parameterization is a special case of the affine Grassmannian
+// manifold (see https://en.wikipedia.org/wiki/Affine_Grassmannian_(manifold))
+// for the case Graff_1(R^n).
+template <int AmbientSpaceDimension>
+class LineParameterization : public LocalParameterization {
+ public:
+ static_assert(AmbientSpaceDimension >= 2,
+ "The ambient space must be at least 2");
+
+ bool Plus(const double* x,
+ const double* delta,
+ double* x_plus_delta) const override;
+ bool ComputeJacobian(const double* x, double* jacobian) const override;
+ int GlobalSize() const override { return 2 * AmbientSpaceDimension; }
+ int LocalSize() const override { return 2 * (AmbientSpaceDimension - 1); }
+};
+
// Construct a local parameterization by taking the Cartesian product
// of a number of other local parameterizations. This is useful, when
// a parameter block is the cartesian product of two or more
// manifolds. For example the parameters of a camera consist of a
// rotation and a translation, i.e., SO(3) x R^3.
//
-// Currently this class supports taking the cartesian product of up to
-// four local parameterizations.
-//
// Example usage:
//
// ProductParameterization product_param(new QuaterionionParameterization(),
@@ -281,35 +304,51 @@
// rotation is represented using a quaternion.
class CERES_EXPORT ProductParameterization : public LocalParameterization {
public:
+ ProductParameterization(const ProductParameterization&) = delete;
+ ProductParameterization& operator=(const ProductParameterization&) = delete;
+ virtual ~ProductParameterization() {}
//
- // NOTE: All the constructors take ownership of the input local
+ // NOTE: The constructor takes ownership of the input local
// parameterizations.
//
- ProductParameterization(LocalParameterization* local_param1,
- LocalParameterization* local_param2);
+ template <typename... LocalParams>
+ ProductParameterization(LocalParams*... local_params)
+ : local_params_(sizeof...(LocalParams)),
+ local_size_{0},
+ global_size_{0},
+ buffer_size_{0} {
+ constexpr int kNumLocalParams = sizeof...(LocalParams);
+ static_assert(kNumLocalParams >= 2,
+ "At least two local parameterizations must be specified.");
- ProductParameterization(LocalParameterization* local_param1,
- LocalParameterization* local_param2,
- LocalParameterization* local_param3);
+ using LocalParameterizationPtr = std::unique_ptr<LocalParameterization>;
- ProductParameterization(LocalParameterization* local_param1,
- LocalParameterization* local_param2,
- LocalParameterization* local_param3,
- LocalParameterization* local_param4);
+ // Wrap all raw pointers into std::unique_ptr for exception safety.
+ std::array<LocalParameterizationPtr, kNumLocalParams> local_params_array{
+ LocalParameterizationPtr(local_params)...};
- virtual ~ProductParameterization();
- virtual bool Plus(const double* x,
- const double* delta,
- double* x_plus_delta) const;
- virtual bool ComputeJacobian(const double* x,
- double* jacobian) const;
- virtual int GlobalSize() const { return global_size_; }
- virtual int LocalSize() const { return local_size_; }
+ // Initialize internal state.
+ for (int i = 0; i < kNumLocalParams; ++i) {
+ LocalParameterizationPtr& param = local_params_[i];
+ param = std::move(local_params_array[i]);
+
+ buffer_size_ =
+ std::max(buffer_size_, param->LocalSize() * param->GlobalSize());
+ global_size_ += param->GlobalSize();
+ local_size_ += param->LocalSize();
+ }
+ }
+
+ bool Plus(const double* x,
+ const double* delta,
+ double* x_plus_delta) const override;
+ bool ComputeJacobian(const double* x,
+ double* jacobian) const override;
+ int GlobalSize() const override { return global_size_; }
+ int LocalSize() const override { return local_size_; }
private:
- void Init();
-
- std::vector<LocalParameterization*> local_params_;
+ std::vector<std::unique_ptr<LocalParameterization>> local_params_;
int local_size_;
int global_size_;
int buffer_size_;
@@ -317,6 +356,8 @@
} // namespace ceres
+// clang-format off
#include "ceres/internal/reenable_warnings.h"
+#include "ceres/internal/line_parameterization.h"
#endif // CERES_PUBLIC_LOCAL_PARAMETERIZATION_H_
diff --git a/include/ceres/loss_function.h b/include/ceres/loss_function.h
index 97a70b6..7aabf7d 100644
--- a/include/ceres/loss_function.h
+++ b/include/ceres/loss_function.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -76,9 +76,10 @@
#define CERES_PUBLIC_LOSS_FUNCTION_H_
#include <memory>
-#include "glog/logging.h"
-#include "ceres/types.h"
+
#include "ceres/internal/disable_warnings.h"
+#include "ceres/types.h"
+#include "glog/logging.h"
namespace ceres {
@@ -118,7 +119,6 @@
// Note: in the region of interest (i.e. s < 3) we have:
// TrivialLoss >= HuberLoss >= SoftLOneLoss >= CauchyLoss
-
// This corresponds to no robustification.
//
// rho(s) = s
@@ -130,7 +130,7 @@
// thing.
class CERES_EXPORT TrivialLoss : public LossFunction {
public:
- virtual void Evaluate(double, double*) const;
+ void Evaluate(double, double*) const override;
};
// Scaling
@@ -173,8 +173,8 @@
// http://en.wikipedia.org/wiki/Huber_Loss_Function
class CERES_EXPORT HuberLoss : public LossFunction {
public:
- explicit HuberLoss(double a) : a_(a), b_(a * a) { }
- virtual void Evaluate(double, double*) const;
+ explicit HuberLoss(double a) : a_(a), b_(a * a) {}
+ void Evaluate(double, double*) const override;
private:
const double a_;
@@ -186,11 +186,11 @@
//
// rho(s) = 2 (sqrt(1 + s) - 1).
//
-// At s = 0: rho = [0, 1, -1/2].
+// At s = 0: rho = [0, 1, -1 / (2 * a^2)].
class CERES_EXPORT SoftLOneLoss : public LossFunction {
public:
- explicit SoftLOneLoss(double a) : b_(a * a), c_(1 / b_) { }
- virtual void Evaluate(double, double*) const;
+ explicit SoftLOneLoss(double a) : b_(a * a), c_(1 / b_) {}
+ void Evaluate(double, double*) const override;
private:
// b = a^2.
@@ -203,11 +203,11 @@
//
// rho(s) = log(1 + s).
//
-// At s = 0: rho = [0, 1, -1].
+// At s = 0: rho = [0, 1, -1 / a^2].
class CERES_EXPORT CauchyLoss : public LossFunction {
public:
- explicit CauchyLoss(double a) : b_(a * a), c_(1 / b_) { }
- virtual void Evaluate(double, double*) const;
+ explicit CauchyLoss(double a) : b_(a * a), c_(1 / b_) {}
+ void Evaluate(double, double*) const override;
private:
// b = a^2.
@@ -227,8 +227,8 @@
// At s = 0: rho = [0, 1, 0].
class CERES_EXPORT ArctanLoss : public LossFunction {
public:
- explicit ArctanLoss(double a) : a_(a), b_(1 / (a * a)) { }
- virtual void Evaluate(double, double*) const;
+ explicit ArctanLoss(double a) : a_(a), b_(1 / (a * a)) {}
+ void Evaluate(double, double*) const override;
private:
const double a_;
@@ -267,7 +267,7 @@
class CERES_EXPORT TolerantLoss : public LossFunction {
public:
explicit TolerantLoss(double a, double b);
- virtual void Evaluate(double, double*) const;
+ void Evaluate(double, double*) const override;
private:
const double a_, b_, c_;
@@ -276,16 +276,17 @@
// This is the Tukey biweight loss function which aggressively
// attempts to suppress large errors.
//
-// The term is computed as:
+// The term is computed as follows where the equations are scaled by a
+// factor of 2 because the cost function is given by 1/2 rho(s):
//
-// rho(s) = a^2 / 6 * (1 - (1 - s / a^2)^3 ) for s <= a^2,
-// rho(s) = a^2 / 6 for s > a^2.
+// rho(s) = a^2 / 3 * (1 - (1 - s / a^2)^3 ) for s <= a^2,
+// rho(s) = a^2 / 3 for s > a^2.
//
-// At s = 0: rho = [0, 0.5, -1 / a^2]
+// At s = 0: rho = [0, 1, -2 / a^2]
class CERES_EXPORT TukeyLoss : public ceres::LossFunction {
public:
- explicit TukeyLoss(double a) : a_squared_(a * a) { }
- virtual void Evaluate(double, double*) const;
+ explicit TukeyLoss(double a) : a_squared_(a * a) {}
+ void Evaluate(double, double*) const override;
private:
const double a_squared_;
@@ -296,10 +297,12 @@
// The loss functions must not be NULL.
class CERES_EXPORT ComposedLoss : public LossFunction {
public:
- explicit ComposedLoss(const LossFunction* f, Ownership ownership_f,
- const LossFunction* g, Ownership ownership_g);
+ explicit ComposedLoss(const LossFunction* f,
+ Ownership ownership_f,
+ const LossFunction* g,
+ Ownership ownership_g);
virtual ~ComposedLoss();
- virtual void Evaluate(double, double*) const;
+ void Evaluate(double, double*) const override;
private:
std::unique_ptr<const LossFunction> f_, g_;
@@ -328,8 +331,8 @@
// Constructs a ScaledLoss wrapping another loss function. Takes
// ownership of the wrapped loss function or not depending on the
// ownership parameter.
- ScaledLoss(const LossFunction* rho, double a, Ownership ownership) :
- rho_(rho), a_(a), ownership_(ownership) { }
+ ScaledLoss(const LossFunction* rho, double a, Ownership ownership)
+ : rho_(rho), a_(a), ownership_(ownership) {}
ScaledLoss(const ScaledLoss&) = delete;
void operator=(const ScaledLoss&) = delete;
@@ -338,7 +341,7 @@
rho_.release();
}
}
- virtual void Evaluate(double, double*) const;
+ void Evaluate(double, double*) const override;
private:
std::unique_ptr<const LossFunction> rho_;
@@ -387,8 +390,7 @@
class CERES_EXPORT LossFunctionWrapper : public LossFunction {
public:
LossFunctionWrapper(LossFunction* rho, Ownership ownership)
- : rho_(rho), ownership_(ownership) {
- }
+ : rho_(rho), ownership_(ownership) {}
LossFunctionWrapper(const LossFunctionWrapper&) = delete;
void operator=(const LossFunctionWrapper&) = delete;
@@ -399,13 +401,12 @@
}
}
- virtual void Evaluate(double sq_norm, double out[3]) const {
+ void Evaluate(double sq_norm, double out[3]) const override {
if (rho_.get() == NULL) {
out[0] = sq_norm;
out[1] = 1.0;
out[2] = 0.0;
- }
- else {
+ } else {
rho_->Evaluate(sq_norm, out);
}
}
diff --git a/include/ceres/normal_prior.h b/include/ceres/normal_prior.h
index cd98b4c..14ab379 100644
--- a/include/ceres/normal_prior.h
+++ b/include/ceres/normal_prior.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -35,8 +35,8 @@
#define CERES_PUBLIC_NORMAL_PRIOR_H_
#include "ceres/cost_function.h"
-#include "ceres/internal/eigen.h"
#include "ceres/internal/disable_warnings.h"
+#include "ceres/internal/eigen.h"
namespace ceres {
@@ -57,15 +57,15 @@
// which would be the case if the covariance matrix S is rank
// deficient.
-class CERES_EXPORT NormalPrior: public CostFunction {
+class CERES_EXPORT NormalPrior : public CostFunction {
public:
// Check that the number of rows in the vector b are the same as the
// number of columns in the matrix A, crash otherwise.
NormalPrior(const Matrix& A, const Vector& b);
+ bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) const override;
- virtual bool Evaluate(double const* const* parameters,
- double* residuals,
- double** jacobians) const;
private:
Matrix A_;
Vector b_;
diff --git a/include/ceres/numeric_diff_cost_function.h b/include/ceres/numeric_diff_cost_function.h
index f6cd32a..cf7971c 100644
--- a/include/ceres/numeric_diff_cost_function.h
+++ b/include/ceres/numeric_diff_cost_function.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -98,6 +98,8 @@
// NumericDiffCostFunction also supports cost functions with a
// runtime-determined number of residuals. For example:
//
+// clang-format off
+//
// CostFunction* cost_function
// = new NumericDiffCostFunction<MyScalarCostFunctor, CENTRAL, DYNAMIC, 2, 2>(
// new CostFunctorWithDynamicNumResiduals(1.0), ^ ^ ^
@@ -109,6 +111,8 @@
// Indicate dynamic number of residuals --------------------+ | |
// Dimension of x ------------------------------------------------+ |
// Dimension of y ---------------------------------------------------+
+// clang-format on
+//
//
// The central difference method is considerably more accurate at the cost of
// twice as many function evaluations than forward difference. Consider using
@@ -182,37 +186,36 @@
Ownership ownership = TAKE_OWNERSHIP,
int num_residuals = kNumResiduals,
const NumericDiffOptions& options = NumericDiffOptions())
- : functor_(functor),
- ownership_(ownership),
- options_(options) {
+ : functor_(functor), ownership_(ownership), options_(options) {
if (kNumResiduals == DYNAMIC) {
SizedCostFunction<kNumResiduals, Ns...>::set_num_residuals(num_residuals);
}
}
- ~NumericDiffCostFunction() {
+ explicit NumericDiffCostFunction(NumericDiffCostFunction&& other)
+ : functor_(std::move(other.functor_)), ownership_(other.ownership_) {}
+
+ virtual ~NumericDiffCostFunction() {
if (ownership_ != TAKE_OWNERSHIP) {
functor_.release();
}
}
- virtual bool Evaluate(double const* const* parameters,
- double* residuals,
- double** jacobians) const {
+ bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) const override {
using internal::FixedArray;
using internal::NumericDiff;
using ParameterDims =
typename SizedCostFunction<kNumResiduals, Ns...>::ParameterDims;
- using Parameters = typename ParameterDims::Parameters;
constexpr int kNumParameters = ParameterDims::kNumParameters;
constexpr int kNumParameterBlocks = ParameterDims::kNumParameterBlocks;
// Get the function value (residuals) at the the point to evaluate.
- if (!internal::VariadicEvaluate<ParameterDims>(*functor_,
- parameters,
- residuals)) {
+ if (!internal::VariadicEvaluate<ParameterDims>(
+ *functor_, parameters, residuals)) {
return false;
}
@@ -223,21 +226,22 @@
// Create a copy of the parameters which will get mutated.
FixedArray<double> parameters_copy(kNumParameters);
std::array<double*, kNumParameterBlocks> parameters_reference_copy =
- ParameterDims::GetUnpackedParameters(parameters_copy.get());
+ ParameterDims::GetUnpackedParameters(parameters_copy.data());
for (int block = 0; block < kNumParameterBlocks; ++block) {
- memcpy(parameters_reference_copy[block], parameters[block],
+ memcpy(parameters_reference_copy[block],
+ parameters[block],
sizeof(double) * ParameterDims::GetDim(block));
}
- internal::EvaluateJacobianForParameterBlocks<ParameterDims>::template Apply<
- method, kNumResiduals>(
- functor_.get(),
- residuals,
- options_,
- SizedCostFunction<kNumResiduals, Ns...>::num_residuals(),
- parameters_reference_copy.data(),
- jacobians);
+ internal::EvaluateJacobianForParameterBlocks<ParameterDims>::
+ template Apply<method, kNumResiduals>(
+ functor_.get(),
+ residuals,
+ options_,
+ SizedCostFunction<kNumResiduals, Ns...>::num_residuals(),
+ parameters_reference_copy.data(),
+ jacobians);
return true;
}
diff --git a/include/ceres/numeric_diff_options.h b/include/ceres/numeric_diff_options.h
index a77be44..64919ed 100644
--- a/include/ceres/numeric_diff_options.h
+++ b/include/ceres/numeric_diff_options.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
diff --git a/include/ceres/ordered_groups.h b/include/ceres/ordered_groups.h
index 0bff41f..954663c 100644
--- a/include/ceres/ordered_groups.h
+++ b/include/ceres/ordered_groups.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -35,6 +35,7 @@
#include <set>
#include <unordered_map>
#include <vector>
+
#include "ceres/internal/port.h"
#include "glog/logging.h"
@@ -161,17 +162,13 @@
// group for every integer.
int GroupSize(const int group) const {
auto it = group_to_elements_.find(group);
- return (it == group_to_elements_.end()) ? 0 : it->second.size();
+ return (it == group_to_elements_.end()) ? 0 : it->second.size();
}
- int NumElements() const {
- return element_to_group_.size();
- }
+ int NumElements() const { return element_to_group_.size(); }
// Number of groups with one or more elements.
- int NumGroups() const {
- return group_to_elements_.size();
- }
+ int NumGroups() const { return group_to_elements_.size(); }
// The first group with one or more elements. Calling this when
// there are no groups with non-zero elements will result in a
@@ -185,9 +182,7 @@
return group_to_elements_;
}
- const std::map<T, int>& element_to_group() const {
- return element_to_group_;
- }
+ const std::map<T, int>& element_to_group() const { return element_to_group_; }
private:
std::map<int, std::set<T>> group_to_elements_;
diff --git a/include/ceres/problem.h b/include/ceres/problem.h
index 503e0fe..add12ea 100644
--- a/include/ceres/problem.h
+++ b/include/ceres/problem.h
@@ -50,6 +50,7 @@
namespace ceres {
class CostFunction;
+class EvaluationCallback;
class LossFunction;
class LocalParameterization;
class Solver;
@@ -114,8 +115,8 @@
//
// Problem problem;
//
-// problem.AddResidualBlock(new MyUnaryCostFunction(...), NULL, x1);
-// problem.AddResidualBlock(new MyBinaryCostFunction(...), NULL, x2, x3);
+// problem.AddResidualBlock(new MyUnaryCostFunction(...), nullptr, x1);
+// problem.AddResidualBlock(new MyBinaryCostFunction(...), nullptr, x2, x3);
//
// Please see cost_function.h for details of the CostFunction object.
class CERES_EXPORT Problem {
@@ -161,18 +162,38 @@
// A Ceres global context to use for solving this problem. This may help to
// reduce computation time as Ceres can reuse expensive objects to create.
- // The context object can be NULL, in which case Ceres may create one.
+ // The context object can be nullptr, in which case Ceres may create one.
//
// Ceres does NOT take ownership of the pointer.
Context* context = nullptr;
+
+ // Using this callback interface, Ceres can notify you when it is
+ // about to evaluate the residuals or jacobians. With the
+ // callback, you can share computation between residual blocks by
+ // doing the shared computation in
+ // EvaluationCallback::PrepareForEvaluation() before Ceres calls
+ // CostFunction::Evaluate(). It also enables caching results
+ // between a pure residual evaluation and a residual & jacobian
+ // evaluation.
+ //
+ // Problem DOES NOT take ownership of the callback.
+ //
+ // NOTE: Evaluation callbacks are incompatible with inner
+ // iterations. So calling Solve with
+ // Solver::Options::use_inner_iterations = true on a Problem with
+ // a non-null evaluation callback is an error.
+ EvaluationCallback* evaluation_callback = nullptr;
};
// The default constructor is equivalent to the
// invocation Problem(Problem::Options()).
Problem();
explicit Problem(const Options& options);
+ Problem(Problem&&);
+ Problem& operator=(Problem&&);
+
Problem(const Problem&) = delete;
- void operator=(const Problem&) = delete;
+ Problem& operator=(const Problem&) = delete;
~Problem();
@@ -181,7 +202,7 @@
// parameter blocks it expects. The function checks that these match
// the sizes of the parameter blocks listed in parameter_blocks. The
// program aborts if a mismatch is detected. loss_function can be
- // NULL, in which case the cost of the term is just the squared norm
+ // nullptr, in which case the cost of the term is just the squared norm
// of the residuals.
//
// The user has the option of explicitly adding the parameter blocks
@@ -210,8 +231,8 @@
//
// Problem problem;
//
- // problem.AddResidualBlock(new MyUnaryCostFunction(...), NULL, x1);
- // problem.AddResidualBlock(new MyBinaryCostFunction(...), NULL, x2, x1);
+ // problem.AddResidualBlock(new MyUnaryCostFunction(...), nullptr, x1);
+ // problem.AddResidualBlock(new MyBinaryCostFunction(...), nullptr, x2, x1);
//
// Add a residual block by listing the parameter block pointers
// directly instead of wapping them in a container.
@@ -221,7 +242,8 @@
double* x0,
Ts*... xs) {
const std::array<double*, sizeof...(Ts) + 1> parameter_blocks{{x0, xs...}};
- return AddResidualBlock(cost_function, loss_function,
+ return AddResidualBlock(cost_function,
+ loss_function,
parameter_blocks.data(),
static_cast<int>(parameter_blocks.size()));
}
@@ -234,11 +256,10 @@
// Add a residual block by providing a pointer to the parameter block array
// and the number of parameter blocks.
- ResidualBlockId AddResidualBlock(
- CostFunction* cost_function,
- LossFunction* loss_function,
- double* const* const parameter_blocks,
- int num_parameter_blocks);
+ ResidualBlockId AddResidualBlock(CostFunction* cost_function,
+ LossFunction* loss_function,
+ double* const* const parameter_blocks,
+ int num_parameter_blocks);
// Add a parameter block with appropriate size to the problem.
// Repeated calls with the same arguments are ignored. Repeated
@@ -268,7 +289,7 @@
// ordering, rendering the jacobian or residuals returned from the solver
// uninterpretable. If you depend on the evaluated jacobian, do not use
// remove! This may change in a future release.
- void RemoveParameterBlock(double* values);
+ void RemoveParameterBlock(const double* values);
// Remove a residual block from the problem. Any parameters that the residual
// block depends on are not removed. The cost and loss functions for the
@@ -282,27 +303,31 @@
void RemoveResidualBlock(ResidualBlockId residual_block);
// Hold the indicated parameter block constant during optimization.
- void SetParameterBlockConstant(double* values);
+ void SetParameterBlockConstant(const double* values);
// Allow the indicated parameter block to vary during optimization.
void SetParameterBlockVariable(double* values);
- // Returns true if a parameter block is set constant, and false otherwise.
- bool IsParameterBlockConstant(double* values) const;
+ // Returns true if a parameter block is set constant, and false
+ // otherwise. A parameter block may be set constant in two ways:
+ // either by calling SetParameterBlockConstant or by associating a
+ // LocalParameterization with a zero dimensional tangent space with
+ // it.
+ bool IsParameterBlockConstant(const double* values) const;
// Set the local parameterization for one of the parameter blocks.
// The local_parameterization is owned by the Problem by default. It
// is acceptable to set the same parameterization for multiple
// parameters; the destructor is careful to delete local
- // parameterizations only once. The local parameterization can only
- // be set once per parameter, and cannot be changed once set.
+ // parameterizations only once. Calling SetParameterization with
+ // nullptr will clear any previously set parameterization.
void SetParameterization(double* values,
LocalParameterization* local_parameterization);
// Get the local parameterization object associated with this
// parameter block. If there is no parameterization object
- // associated then NULL is returned.
- const LocalParameterization* GetParameterization(double* values) const;
+ // associated then nullptr is returned.
+ const LocalParameterization* GetParameterization(const double* values) const;
// Set the lower/upper bound for the parameter at position "index".
void SetParameterLowerBound(double* values, int index, double lower_bound);
@@ -312,8 +337,8 @@
// "index". If the parameter is not bounded by the user, then its
// lower bound is -std::numeric_limits<double>::max() and upper
// bound is std::numeric_limits<double>::max().
- double GetParameterLowerBound(double* values, int index) const;
- double GetParameterUpperBound(double* values, int index) const;
+ double GetParameterLowerBound(const double* values, int index) const;
+ double GetParameterUpperBound(const double* values, int index) const;
// Number of parameter blocks in the problem. Always equals
// parameter_blocks().size() and parameter_block_sizes().size().
@@ -361,7 +386,7 @@
const CostFunction* GetCostFunctionForResidualBlock(
const ResidualBlockId residual_block) const;
- // Get the LossFunction for the given residual block. Returns NULL
+ // Get the LossFunction for the given residual block. Returns nullptr
// if no loss function is associated with this residual block.
const LossFunction* GetLossFunctionForResidualBlock(
const ResidualBlockId residual_block) const;
@@ -415,7 +440,7 @@
int num_threads = 1;
};
- // Evaluate Problem. Any of the output pointers can be NULL. Which
+ // Evaluate Problem. Any of the output pointers can be nullptr. Which
// residual blocks and parameter blocks are used is controlled by
// the EvaluateOptions struct above.
//
@@ -425,16 +450,18 @@
//
// Problem problem;
// double x = 1;
- // problem.AddResidualBlock(new MyCostFunction, NULL, &x);
+ // problem.AddResidualBlock(new MyCostFunction, nullptr, &x);
//
// double cost = 0.0;
- // problem.Evaluate(Problem::EvaluateOptions(), &cost, NULL, NULL, NULL);
+ // problem.Evaluate(Problem::EvaluateOptions(), &cost,
+ // nullptr, nullptr, nullptr);
//
// The cost is evaluated at x = 1. If you wish to evaluate the
// problem at x = 2, then
//
- // x = 2;
- // problem.Evaluate(Problem::EvaluateOptions(), &cost, NULL, NULL, NULL);
+ // x = 2;
+ // problem.Evaluate(Problem::EvaluateOptions(), &cost,
+ // nullptr, nullptr, nullptr);
//
// is the way to do so.
//
@@ -448,16 +475,81 @@
// Note 3: This function cannot be called while the problem is being
// solved, for example it cannot be called from an IterationCallback
// at the end of an iteration during a solve.
+ //
+ // Note 4: If an EvaluationCallback is associated with the problem,
+ // then its PrepareForEvaluation method will be called every time
+ // this method is called with new_point = true.
bool Evaluate(const EvaluateOptions& options,
double* cost,
std::vector<double>* residuals,
std::vector<double>* gradient,
CRSMatrix* jacobian);
+ // Evaluates the residual block, storing the scalar cost in *cost,
+ // the residual components in *residuals, and the jacobians between
+ // the parameters and residuals in jacobians[i], in row-major order.
+ //
+ // If residuals is nullptr, the residuals are not computed.
+ //
+ // If jacobians is nullptr, no Jacobians are computed. If
+ // jacobians[i] is nullptr, then the Jacobian for that parameter
+ // block is not computed.
+ //
+ // It is not okay to request the Jacobian w.r.t a parameter block
+ // that is constant.
+ //
+ // The return value indicates the success or failure. Even if the
+ // function returns false, the caller should expect the output
+ // memory locations to have been modified.
+ //
+ // The returned cost and jacobians have had robustification and
+ // local parameterizations applied already; for example, the
+ // jacobian for a 4-dimensional quaternion parameter using the
+ // "QuaternionParameterization" is num_residuals by 3 instead of
+ // num_residuals by 4.
+ //
+ // apply_loss_function as the name implies allows the user to switch
+ // the application of the loss function on and off.
+ //
+ // If an EvaluationCallback is associated with the problem, then its
+ // PrepareForEvaluation method will be called every time this method
+ // is called with new_point = true. This conservatively assumes that
+ // the user may have changed the parameter values since the previous
+ // call to evaluate / solve. For improved efficiency, and only if
+ // you know that the parameter values have not changed between
+ // calls, see EvaluateResidualBlockAssumingParametersUnchanged().
+ bool EvaluateResidualBlock(ResidualBlockId residual_block_id,
+ bool apply_loss_function,
+ double* cost,
+ double* residuals,
+ double** jacobians) const;
+
+ // Same as EvaluateResidualBlock except that if an
+ // EvaluationCallback is associated with the problem, then its
+ // PrepareForEvaluation method will be called every time this method
+ // is called with new_point = false.
+ //
+ // This means, if an EvaluationCallback is associated with the
+ // problem then it is the user's responsibility to call
+ // PrepareForEvaluation before calling this method if necessary,
+ // i.e. iff the parameter values have been changed since the last
+ // call to evaluate / solve.'
+ //
+ // This is because, as the name implies, we assume that the
+ // parameter blocks did not change since the last time
+ // PrepareForEvaluation was called (via Solve, Evaluate or
+ // EvaluateResidualBlock).
+ bool EvaluateResidualBlockAssumingParametersUnchanged(
+ ResidualBlockId residual_block_id,
+ bool apply_loss_function,
+ double* cost,
+ double* residuals,
+ double** jacobians) const;
+
private:
friend class Solver;
friend class Covariance;
- std::unique_ptr<internal::ProblemImpl> problem_impl_;
+ std::unique_ptr<internal::ProblemImpl> impl_;
};
} // namespace ceres
diff --git a/include/ceres/rotation.h b/include/ceres/rotation.h
index a0530dd..0c82a41 100644
--- a/include/ceres/rotation.h
+++ b/include/ceres/rotation.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -49,6 +49,8 @@
#include <cmath>
#include <limits>
+#include "glog/logging.h"
+
namespace ceres {
// Trivial wrapper to index linear arrays as matrices, given a fixed
@@ -82,7 +84,7 @@
// and quaternion is a 4-tuple that will contain the resulting quaternion.
// The implementation may be used with auto-differentiation up to the first
// derivative, higher derivatives may have unexpected results near the origin.
-template<typename T>
+template <typename T>
void AngleAxisToQuaternion(const T* angle_axis, T* quaternion);
// Convert a quaternion to the equivalent combined axis-angle representation.
@@ -91,7 +93,7 @@
// rotation in radians, and whose direction is the axis of rotation.
// The implementation may be used with auto-differentiation up to the first
// derivative, higher derivatives may have unexpected results near the origin.
-template<typename T>
+template <typename T>
void QuaternionToAngleAxis(const T* quaternion, T* angle_axis);
// Conversions between 3x3 rotation matrix (in column major order) and
@@ -102,8 +104,7 @@
template <typename T, int row_stride, int col_stride>
void RotationMatrixToQuaternion(
- const MatrixAdapter<const T, row_stride, col_stride>& R,
- T* quaternion);
+ const MatrixAdapter<const T, row_stride, col_stride>& R, T* quaternion);
// Conversions between 3x3 rotation matrix (in column major order) and
// axis-angle rotation representations. Templated for use with
@@ -113,16 +114,14 @@
template <typename T, int row_stride, int col_stride>
void RotationMatrixToAngleAxis(
- const MatrixAdapter<const T, row_stride, col_stride>& R,
- T* angle_axis);
+ const MatrixAdapter<const T, row_stride, col_stride>& R, T* angle_axis);
template <typename T>
void AngleAxisToRotationMatrix(const T* angle_axis, T* R);
template <typename T, int row_stride, int col_stride>
void AngleAxisToRotationMatrix(
- const T* angle_axis,
- const MatrixAdapter<T, row_stride, col_stride>& R);
+ const T* angle_axis, const MatrixAdapter<T, row_stride, col_stride>& R);
// Conversions between 3x3 rotation matrix (in row major order) and
// Euler angle (in degrees) rotation representations.
@@ -135,8 +134,7 @@
template <typename T, int row_stride, int col_stride>
void EulerAnglesToRotationMatrix(
- const T* euler,
- const MatrixAdapter<T, row_stride, col_stride>& R);
+ const T* euler, const MatrixAdapter<T, row_stride, col_stride>& R);
// Convert a 4-vector to a 3x3 scaled rotation matrix.
//
@@ -157,25 +155,23 @@
// such that det(Q) = 1 and Q*Q' = I
//
// WARNING: The rotation matrix is ROW MAJOR
-template <typename T> inline
-void QuaternionToScaledRotation(const T q[4], T R[3 * 3]);
+template <typename T>
+inline void QuaternionToScaledRotation(const T q[4], T R[3 * 3]);
-template <typename T, int row_stride, int col_stride> inline
-void QuaternionToScaledRotation(
- const T q[4],
- const MatrixAdapter<T, row_stride, col_stride>& R);
+template <typename T, int row_stride, int col_stride>
+inline void QuaternionToScaledRotation(
+ const T q[4], const MatrixAdapter<T, row_stride, col_stride>& R);
// Same as above except that the rotation matrix is normalized by the
// Frobenius norm, so that R * R' = I (and det(R) = 1).
//
// WARNING: The rotation matrix is ROW MAJOR
-template <typename T> inline
-void QuaternionToRotation(const T q[4], T R[3 * 3]);
+template <typename T>
+inline void QuaternionToRotation(const T q[4], T R[3 * 3]);
-template <typename T, int row_stride, int col_stride> inline
-void QuaternionToRotation(
- const T q[4],
- const MatrixAdapter<T, row_stride, col_stride>& R);
+template <typename T, int row_stride, int col_stride>
+inline void QuaternionToRotation(
+ const T q[4], const MatrixAdapter<T, row_stride, col_stride>& R);
// Rotates a point pt by a quaternion q:
//
@@ -185,37 +181,54 @@
// write the transform as (something)*pt + pt, as is clear from the
// formula below. If you pass in a quaternion with |q|^2 = 2 then you
// WILL NOT get back 2 times the result you get for a unit quaternion.
-template <typename T> inline
-void UnitQuaternionRotatePoint(const T q[4], const T pt[3], T result[3]);
+//
+// Inplace rotation is not supported. pt and result must point to different
+// memory locations, otherwise the result will be undefined.
+template <typename T>
+inline void UnitQuaternionRotatePoint(const T q[4], const T pt[3], T result[3]);
// With this function you do not need to assume that q has unit norm.
// It does assume that the norm is non-zero.
-template <typename T> inline
-void QuaternionRotatePoint(const T q[4], const T pt[3], T result[3]);
+//
+// Inplace rotation is not supported. pt and result must point to different
+// memory locations, otherwise the result will be undefined.
+template <typename T>
+inline void QuaternionRotatePoint(const T q[4], const T pt[3], T result[3]);
// zw = z * w, where * is the Quaternion product between 4 vectors.
-template<typename T> inline
-void QuaternionProduct(const T z[4], const T w[4], T zw[4]);
+//
+// Inplace quaternion product is not supported. The resulting quaternion zw must
+// not share the memory with the input quaternion z and w, otherwise the result
+// will be undefined.
+template <typename T>
+inline void QuaternionProduct(const T z[4], const T w[4], T zw[4]);
// xy = x cross y;
-template<typename T> inline
-void CrossProduct(const T x[3], const T y[3], T x_cross_y[3]);
+//
+// Inplace cross product is not supported. The resulting vector x_cross_y must
+// not share the memory with the input vectors x and y, otherwise the result
+// will be undefined.
+template <typename T>
+inline void CrossProduct(const T x[3], const T y[3], T x_cross_y[3]);
-template<typename T> inline
-T DotProduct(const T x[3], const T y[3]);
+template <typename T>
+inline T DotProduct(const T x[3], const T y[3]);
// y = R(angle_axis) * x;
-template<typename T> inline
-void AngleAxisRotatePoint(const T angle_axis[3], const T pt[3], T result[3]);
+//
+// Inplace rotation is not supported. pt and result must point to different
+// memory locations, otherwise the result will be undefined.
+template <typename T>
+inline void AngleAxisRotatePoint(const T angle_axis[3],
+ const T pt[3],
+ T result[3]);
// --- IMPLEMENTATION
-template<typename T, int row_stride, int col_stride>
+template <typename T, int row_stride, int col_stride>
struct MatrixAdapter {
T* pointer_;
- explicit MatrixAdapter(T* pointer)
- : pointer_(pointer)
- {}
+ explicit MatrixAdapter(T* pointer) : pointer_(pointer) {}
T& operator()(int r, int c) const {
return pointer_[r * row_stride + c * col_stride];
@@ -232,7 +245,7 @@
return MatrixAdapter<T, 3, 1>(pointer);
}
-template<typename T>
+template <typename T>
inline void AngleAxisToQuaternion(const T* angle_axis, T* quaternion) {
const T& a0 = angle_axis[0];
const T& a1 = angle_axis[1];
@@ -261,7 +274,7 @@
}
}
-template<typename T>
+template <typename T>
inline void QuaternionToAngleAxis(const T* quaternion, T* angle_axis) {
const T& q1 = quaternion[1];
const T& q2 = quaternion[2];
@@ -288,9 +301,8 @@
// = atan(-sin(theta), -cos(theta))
//
const T two_theta =
- T(2.0) * ((cos_theta < 0.0)
- ? atan2(-sin_theta, -cos_theta)
- : atan2(sin_theta, cos_theta));
+ T(2.0) * ((cos_theta < T(0.0)) ? atan2(-sin_theta, -cos_theta)
+ : atan2(sin_theta, cos_theta));
const T k = two_theta / sin_theta;
angle_axis[0] = q1 * k;
angle_axis[1] = q2 * k;
@@ -308,16 +320,15 @@
}
template <typename T>
-void RotationMatrixToQuaternion(const T* R, T* angle_axis) {
- RotationMatrixToQuaternion(ColumnMajorAdapter3x3(R), angle_axis);
+void RotationMatrixToQuaternion(const T* R, T* quaternion) {
+ RotationMatrixToQuaternion(ColumnMajorAdapter3x3(R), quaternion);
}
// This algorithm comes from "Quaternion Calculus and Fast Animation",
// Ken Shoemake, 1987 SIGGRAPH course notes
template <typename T, int row_stride, int col_stride>
void RotationMatrixToQuaternion(
- const MatrixAdapter<const T, row_stride, col_stride>& R,
- T* quaternion) {
+ const MatrixAdapter<const T, row_stride, col_stride>& R, T* quaternion) {
const T trace = R(0, 0) + R(1, 1) + R(2, 2);
if (trace >= 0.0) {
T t = sqrt(trace + T(1.0));
@@ -359,8 +370,7 @@
template <typename T, int row_stride, int col_stride>
void RotationMatrixToAngleAxis(
- const MatrixAdapter<const T, row_stride, col_stride>& R,
- T* angle_axis) {
+ const MatrixAdapter<const T, row_stride, col_stride>& R, T* angle_axis) {
T quaternion[4];
RotationMatrixToQuaternion(R, quaternion);
QuaternionToAngleAxis(quaternion, angle_axis);
@@ -374,8 +384,7 @@
template <typename T, int row_stride, int col_stride>
void AngleAxisToRotationMatrix(
- const T* angle_axis,
- const MatrixAdapter<T, row_stride, col_stride>& R) {
+ const T* angle_axis, const MatrixAdapter<T, row_stride, col_stride>& R) {
static const T kOne = T(1.0);
const T theta2 = DotProduct(angle_axis, angle_axis);
if (theta2 > T(std::numeric_limits<double>::epsilon())) {
@@ -390,6 +399,7 @@
const T costheta = cos(theta);
const T sintheta = sin(theta);
+ // clang-format off
R(0, 0) = costheta + wx*wx*(kOne - costheta);
R(1, 0) = wz*sintheta + wx*wy*(kOne - costheta);
R(2, 0) = -wy*sintheta + wx*wz*(kOne - costheta);
@@ -399,15 +409,16 @@
R(0, 2) = wy*sintheta + wx*wz*(kOne - costheta);
R(1, 2) = -wx*sintheta + wy*wz*(kOne - costheta);
R(2, 2) = costheta + wz*wz*(kOne - costheta);
+ // clang-format on
} else {
// Near zero, we switch to using the first order Taylor expansion.
- R(0, 0) = kOne;
- R(1, 0) = angle_axis[2];
+ R(0, 0) = kOne;
+ R(1, 0) = angle_axis[2];
R(2, 0) = -angle_axis[1];
R(0, 1) = -angle_axis[2];
- R(1, 1) = kOne;
- R(2, 1) = angle_axis[0];
- R(0, 2) = angle_axis[1];
+ R(1, 1) = kOne;
+ R(2, 1) = angle_axis[0];
+ R(0, 2) = angle_axis[1];
R(1, 2) = -angle_axis[0];
R(2, 2) = kOne;
}
@@ -422,8 +433,7 @@
template <typename T, int row_stride, int col_stride>
void EulerAnglesToRotationMatrix(
- const T* euler,
- const MatrixAdapter<T, row_stride, col_stride>& R) {
+ const T* euler, const MatrixAdapter<T, row_stride, col_stride>& R) {
const double kPi = 3.14159265358979323846;
const T degrees_to_radians(kPi / 180.0);
@@ -438,28 +448,27 @@
const T c3 = cos(pitch);
const T s3 = sin(pitch);
- R(0, 0) = c1*c2;
- R(0, 1) = -s1*c3 + c1*s2*s3;
- R(0, 2) = s1*s3 + c1*s2*c3;
+ R(0, 0) = c1 * c2;
+ R(0, 1) = -s1 * c3 + c1 * s2 * s3;
+ R(0, 2) = s1 * s3 + c1 * s2 * c3;
- R(1, 0) = s1*c2;
- R(1, 1) = c1*c3 + s1*s2*s3;
- R(1, 2) = -c1*s3 + s1*s2*c3;
+ R(1, 0) = s1 * c2;
+ R(1, 1) = c1 * c3 + s1 * s2 * s3;
+ R(1, 2) = -c1 * s3 + s1 * s2 * c3;
R(2, 0) = -s2;
- R(2, 1) = c2*s3;
- R(2, 2) = c2*c3;
+ R(2, 1) = c2 * s3;
+ R(2, 2) = c2 * c3;
}
-template <typename T> inline
-void QuaternionToScaledRotation(const T q[4], T R[3 * 3]) {
+template <typename T>
+inline void QuaternionToScaledRotation(const T q[4], T R[3 * 3]) {
QuaternionToScaledRotation(q, RowMajorAdapter3x3(R));
}
-template <typename T, int row_stride, int col_stride> inline
-void QuaternionToScaledRotation(
- const T q[4],
- const MatrixAdapter<T, row_stride, col_stride>& R) {
+template <typename T, int row_stride, int col_stride>
+inline void QuaternionToScaledRotation(
+ const T q[4], const MatrixAdapter<T, row_stride, col_stride>& R) {
// Make convenient names for elements of q.
T a = q[0];
T b = q[1];
@@ -478,22 +487,24 @@
T cd = c * d;
T dd = d * d;
- R(0, 0) = aa + bb - cc - dd; R(0, 1) = T(2) * (bc - ad); R(0, 2) = T(2) * (ac + bd); // NOLINT
- R(1, 0) = T(2) * (ad + bc); R(1, 1) = aa - bb + cc - dd; R(1, 2) = T(2) * (cd - ab); // NOLINT
- R(2, 0) = T(2) * (bd - ac); R(2, 1) = T(2) * (ab + cd); R(2, 2) = aa - bb - cc + dd; // NOLINT
+ // clang-format off
+ R(0, 0) = aa + bb - cc - dd; R(0, 1) = T(2) * (bc - ad); R(0, 2) = T(2) * (ac + bd);
+ R(1, 0) = T(2) * (ad + bc); R(1, 1) = aa - bb + cc - dd; R(1, 2) = T(2) * (cd - ab);
+ R(2, 0) = T(2) * (bd - ac); R(2, 1) = T(2) * (ab + cd); R(2, 2) = aa - bb - cc + dd;
+ // clang-format on
}
-template <typename T> inline
-void QuaternionToRotation(const T q[4], T R[3 * 3]) {
+template <typename T>
+inline void QuaternionToRotation(const T q[4], T R[3 * 3]) {
QuaternionToRotation(q, RowMajorAdapter3x3(R));
}
-template <typename T, int row_stride, int col_stride> inline
-void QuaternionToRotation(const T q[4],
- const MatrixAdapter<T, row_stride, col_stride>& R) {
+template <typename T, int row_stride, int col_stride>
+inline void QuaternionToRotation(
+ const T q[4], const MatrixAdapter<T, row_stride, col_stride>& R) {
QuaternionToScaledRotation(q, R);
- T normalizer = q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3];
+ T normalizer = q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3];
normalizer = T(1) / normalizer;
for (int i = 0; i < 3; ++i) {
@@ -503,8 +514,13 @@
}
}
-template <typename T> inline
-void UnitQuaternionRotatePoint(const T q[4], const T pt[3], T result[3]) {
+template <typename T>
+inline void UnitQuaternionRotatePoint(const T q[4],
+ const T pt[3],
+ T result[3]) {
+ DCHECK_NE(pt, result) << "Inplace rotation is not supported.";
+
+ // clang-format off
const T t2 = q[0] * q[1];
const T t3 = q[0] * q[2];
const T t4 = q[0] * q[3];
@@ -517,50 +533,63 @@
result[0] = T(2) * ((t8 + t1) * pt[0] + (t6 - t4) * pt[1] + (t3 + t7) * pt[2]) + pt[0]; // NOLINT
result[1] = T(2) * ((t4 + t6) * pt[0] + (t5 + t1) * pt[1] + (t9 - t2) * pt[2]) + pt[1]; // NOLINT
result[2] = T(2) * ((t7 - t3) * pt[0] + (t2 + t9) * pt[1] + (t5 + t8) * pt[2]) + pt[2]; // NOLINT
+ // clang-format on
}
-template <typename T> inline
-void QuaternionRotatePoint(const T q[4], const T pt[3], T result[3]) {
+template <typename T>
+inline void QuaternionRotatePoint(const T q[4], const T pt[3], T result[3]) {
+ DCHECK_NE(pt, result) << "Inplace rotation is not supported.";
+
// 'scale' is 1 / norm(q).
- const T scale = T(1) / sqrt(q[0] * q[0] +
- q[1] * q[1] +
- q[2] * q[2] +
- q[3] * q[3]);
+ const T scale =
+ T(1) / sqrt(q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3]);
// Make unit-norm version of q.
const T unit[4] = {
- scale * q[0],
- scale * q[1],
- scale * q[2],
- scale * q[3],
+ scale * q[0],
+ scale * q[1],
+ scale * q[2],
+ scale * q[3],
};
UnitQuaternionRotatePoint(unit, pt, result);
}
-template<typename T> inline
-void QuaternionProduct(const T z[4], const T w[4], T zw[4]) {
+template <typename T>
+inline void QuaternionProduct(const T z[4], const T w[4], T zw[4]) {
+ DCHECK_NE(z, zw) << "Inplace quaternion product is not supported.";
+ DCHECK_NE(w, zw) << "Inplace quaternion product is not supported.";
+
+ // clang-format off
zw[0] = z[0] * w[0] - z[1] * w[1] - z[2] * w[2] - z[3] * w[3];
zw[1] = z[0] * w[1] + z[1] * w[0] + z[2] * w[3] - z[3] * w[2];
zw[2] = z[0] * w[2] - z[1] * w[3] + z[2] * w[0] + z[3] * w[1];
zw[3] = z[0] * w[3] + z[1] * w[2] - z[2] * w[1] + z[3] * w[0];
+ // clang-format on
}
// xy = x cross y;
-template<typename T> inline
-void CrossProduct(const T x[3], const T y[3], T x_cross_y[3]) {
+template <typename T>
+inline void CrossProduct(const T x[3], const T y[3], T x_cross_y[3]) {
+ DCHECK_NE(x, x_cross_y) << "Inplace cross product is not supported.";
+ DCHECK_NE(y, x_cross_y) << "Inplace cross product is not supported.";
+
x_cross_y[0] = x[1] * y[2] - x[2] * y[1];
x_cross_y[1] = x[2] * y[0] - x[0] * y[2];
x_cross_y[2] = x[0] * y[1] - x[1] * y[0];
}
-template<typename T> inline
-T DotProduct(const T x[3], const T y[3]) {
+template <typename T>
+inline T DotProduct(const T x[3], const T y[3]) {
return (x[0] * y[0] + x[1] * y[1] + x[2] * y[2]);
}
-template<typename T> inline
-void AngleAxisRotatePoint(const T angle_axis[3], const T pt[3], T result[3]) {
+template <typename T>
+inline void AngleAxisRotatePoint(const T angle_axis[3],
+ const T pt[3],
+ T result[3]) {
+ DCHECK_NE(pt, result) << "Inplace rotation is not supported.";
+
const T theta2 = DotProduct(angle_axis, angle_axis);
if (theta2 > T(std::numeric_limits<double>::epsilon())) {
// Away from zero, use the rodriguez formula
@@ -578,15 +607,15 @@
const T sintheta = sin(theta);
const T theta_inverse = T(1.0) / theta;
- const T w[3] = { angle_axis[0] * theta_inverse,
- angle_axis[1] * theta_inverse,
- angle_axis[2] * theta_inverse };
+ const T w[3] = {angle_axis[0] * theta_inverse,
+ angle_axis[1] * theta_inverse,
+ angle_axis[2] * theta_inverse};
// Explicitly inlined evaluation of the cross product for
// performance reasons.
- const T w_cross_pt[3] = { w[1] * pt[2] - w[2] * pt[1],
- w[2] * pt[0] - w[0] * pt[2],
- w[0] * pt[1] - w[1] * pt[0] };
+ const T w_cross_pt[3] = {w[1] * pt[2] - w[2] * pt[1],
+ w[2] * pt[0] - w[0] * pt[2],
+ w[0] * pt[1] - w[1] * pt[0]};
const T tmp =
(w[0] * pt[0] + w[1] * pt[1] + w[2] * pt[2]) * (T(1.0) - costheta);
@@ -611,9 +640,9 @@
//
// Explicitly inlined evaluation of the cross product for
// performance reasons.
- const T w_cross_pt[3] = { angle_axis[1] * pt[2] - angle_axis[2] * pt[1],
- angle_axis[2] * pt[0] - angle_axis[0] * pt[2],
- angle_axis[0] * pt[1] - angle_axis[1] * pt[0] };
+ const T w_cross_pt[3] = {angle_axis[1] * pt[2] - angle_axis[2] * pt[1],
+ angle_axis[2] * pt[0] - angle_axis[0] * pt[2],
+ angle_axis[0] * pt[1] - angle_axis[1] * pt[0]};
result[0] = pt[0] + w_cross_pt[0];
result[1] = pt[1] + w_cross_pt[1];
diff --git a/include/ceres/sized_cost_function.h b/include/ceres/sized_cost_function.h
index 50d0363..8e92f1b 100644
--- a/include/ceres/sized_cost_function.h
+++ b/include/ceres/sized_cost_function.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -53,7 +53,7 @@
static_assert(internal::StaticParameterDims<Ns...>::kIsValid,
"Invalid parameter block dimension detected. Each parameter "
"block dimension must be bigger than zero.");
-
+
using ParameterDims = internal::StaticParameterDims<Ns...>;
SizedCostFunction() {
@@ -61,7 +61,7 @@
*mutable_parameter_block_sizes() = std::vector<int32_t>{Ns...};
}
- virtual ~SizedCostFunction() { }
+ virtual ~SizedCostFunction() {}
// Subclasses must implement Evaluate().
};
diff --git a/include/ceres/solver.h b/include/ceres/solver.h
index 83077e2..61b8dd5 100644
--- a/include/ceres/solver.h
+++ b/include/ceres/solver.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -34,19 +34,19 @@
#include <cmath>
#include <memory>
#include <string>
+#include <unordered_set>
#include <vector>
+
#include "ceres/crs_matrix.h"
-#include "ceres/evaluation_callback.h"
#include "ceres/internal/disable_warnings.h"
#include "ceres/internal/port.h"
#include "ceres/iteration_callback.h"
#include "ceres/ordered_groups.h"
+#include "ceres/problem.h"
#include "ceres/types.h"
namespace ceres {
-class Problem;
-
// Interface for non-linear least squares solvers.
class CERES_EXPORT Solver {
public:
@@ -118,7 +118,7 @@
// method, please see:
//
// Nocedal, J. (1980). "Updating Quasi-Newton Matrices with
- // Limited Storage". Mathematics of Computation 35 (151): 773–782.
+ // Limited Storage". Mathematics of Computation 35 (151): 773-782.
int max_lbfgs_rank = 20;
// As part of the (L)BFGS update step (BFGS) / right-multiply step (L-BFGS),
@@ -188,9 +188,15 @@
//
double min_line_search_step_contraction = 0.6;
- // Maximum number of trial step size iterations during each line search,
- // if a step size satisfying the search conditions cannot be found within
- // this number of trials, the line search will terminate.
+ // Maximum number of trial step size iterations during each line
+ // search, if a step size satisfying the search conditions cannot
+ // be found within this number of trials, the line search will
+ // terminate.
+
+ // The minimum allowed value is 0 for trust region minimizer and 1
+ // otherwise. If 0 is specified for the trust region minimizer,
+ // then line search will not be used when solving constrained
+ // optimization problems.
int max_num_line_search_step_size_iterations = 20;
// Maximum number of restarts of the line search direction algorithm before
@@ -332,6 +338,31 @@
// preconditioner_type is CLUSTER_JACOBI or CLUSTER_TRIDIAGONAL.
VisibilityClusteringType visibility_clustering_type = CANONICAL_VIEWS;
+ // Subset preconditioner is a preconditioner for problems with
+ // general sparsity. Given a subset of residual blocks of a
+ // problem, it uses the corresponding subset of the rows of the
+ // Jacobian to construct a preconditioner.
+ //
+ // Suppose the Jacobian J has been horizontally partitioned as
+ //
+ // J = [P]
+ // [Q]
+ //
+ // Where, Q is the set of rows corresponding to the residual
+ // blocks in residual_blocks_for_subset_preconditioner.
+ //
+ // The preconditioner is the inverse of the matrix Q'Q.
+ //
+ // Obviously, the efficacy of the preconditioner depends on how
+ // well the matrix Q approximates J'J, or how well the chosen
+ // residual blocks approximate the non-linear least squares
+ // problem.
+ //
+ // If Solver::Options::preconditioner_type == SUBSET, then
+ // residual_blocks_for_subset_preconditioner must be non-empty.
+ std::unordered_set<ResidualBlockId>
+ residual_blocks_for_subset_preconditioner;
+
// Ceres supports using multiple dense linear algebra libraries
// for dense matrix factorizations. Currently EIGEN and LAPACK are
// the valid choices. EIGEN is always available, LAPACK refers to
@@ -352,12 +383,12 @@
SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type =
#if !defined(CERES_NO_SUITESPARSE)
SUITE_SPARSE;
-#elif !defined(CERES_NO_ACCELERATE_SPARSE)
- ACCELERATE_SPARSE;
-#elif !defined(CERES_NO_CXSPARSE)
- CX_SPARSE;
#elif defined(CERES_USE_EIGEN_SPARSE)
EIGEN_SPARSE;
+#elif !defined(CERES_NO_CXSPARSE)
+ CX_SPARSE;
+#elif !defined(CERES_NO_ACCELERATE_SPARSE)
+ ACCELERATE_SPARSE;
#else
NO_SPARSE;
#endif
@@ -691,26 +722,18 @@
// optimistic. This number should be exposed for users to change.
double gradient_check_numeric_derivative_relative_step_size = 1e-6;
- // If true, the user's parameter blocks are updated at the end of
- // every Minimizer iteration, otherwise they are updated when the
- // Minimizer terminates. This is useful if, for example, the user
- // wishes to visualize the state of the optimization every iteration
- // (in combination with an IterationCallback).
- //
- // NOTE: If an evaluation_callback is provided, then the behaviour
- // of this flag is slightly different in each case:
- //
- // (1) If update_state_every_iteration = false, then the user's
- // state is changed at every residual and/or jacobian evaluation.
- // Any user provided IterationCallbacks should NOT inspect and
- // depend on the user visible state while the solver is running,
- // since there will be undefined contents.
- //
- // (2) If update_state_every_iteration is true, then the user's
- // state is changed at every residual and/or jacobian evaluation,
- // BUT the solver will ensure that before the user provided
- // IterationCallbacks are called, the user visible state will be
- // updated to the current best point found by the solver.
+ // If update_state_every_iteration is true, then Ceres Solver will
+ // guarantee that at the end of every iteration and before any
+ // user provided IterationCallback is called, the parameter blocks
+ // are updated to the current best solution found by the
+ // solver. Thus the IterationCallback can inspect the values of
+ // the parameter blocks for purposes of computation, visualization
+ // or termination.
+
+ // If update_state_every_iteration is false then there is no such
+ // guarantee, and user provided IterationCallbacks should not
+ // expect to look at the parameter blocks and interpret their
+ // values.
bool update_state_every_iteration = false;
// Callbacks that are executed at the end of each iteration of the
@@ -729,19 +752,6 @@
//
// The solver does NOT take ownership of these pointers.
std::vector<IterationCallback*> callbacks;
-
- // If non-NULL, gets notified when Ceres is about to evaluate the
- // residuals and/or Jacobians. This enables sharing computation
- // between residuals, which in some cases is important for efficient
- // cost evaluation. See evaluation_callback.h for details.
- //
- // NOTE: Evaluation callbacks are incompatible with inner iterations.
- //
- // WARNING: This interacts with update_state_every_iteration. See
- // the documentation for that option for more details.
- //
- // The solver does NOT take ownership of the pointer.
- EvaluationCallback* evaluation_callback = nullptr;
};
struct CERES_EXPORT Summary {
@@ -783,22 +793,22 @@
// accepted. Unless use_non_monotonic_steps is true this is also
// the number of steps in which the objective function value/cost
// went down.
- int num_successful_steps = -1.0;
+ int num_successful_steps = -1;
// Number of minimizer iterations in which the step was rejected
// either because it did not reduce the cost enough or the step
// was not numerically valid.
- int num_unsuccessful_steps = -1.0;
+ int num_unsuccessful_steps = -1;
// Number of times inner iterations were performed.
- int num_inner_iteration_steps = -1.0;
+ int num_inner_iteration_steps = -1;
// Total number of iterations inside the line search algorithm
// across all invocations. We call these iterations "steps" to
// distinguish them from the outer iterations of the line search
// and trust region minimizer algorithms which call the line
// search algorithm as a subroutine.
- int num_line_search_steps = -1.0;
+ int num_line_search_steps = -1;
// All times reported below are wall times.
@@ -829,7 +839,7 @@
int num_linear_solves = -1;
// Time (in seconds) spent evaluating the residual vector.
- double residual_evaluation_time_in_seconds = 1.0;
+ double residual_evaluation_time_in_seconds = -1.0;
// Number of residual only evaluations.
int num_residual_evaluations = -1;
@@ -1043,7 +1053,8 @@
};
// Helper function which avoids going through the interface.
-CERES_EXPORT void Solve(const Solver::Options& options, Problem* problem,
+CERES_EXPORT void Solve(const Solver::Options& options,
+ Problem* problem,
Solver::Summary* summary);
} // namespace ceres
diff --git a/include/ceres/tiny_solver.h b/include/ceres/tiny_solver.h
index b6484fe..47db582 100644
--- a/include/ceres/tiny_solver.h
+++ b/include/ceres/tiny_solver.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2017 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -124,13 +124,17 @@
//
// int NumParameters() const;
//
-template<typename Function,
- typename LinearSolver = Eigen::LDLT<
- Eigen::Matrix<typename Function::Scalar,
- Function::NUM_PARAMETERS,
- Function::NUM_PARAMETERS>>>
+template <typename Function,
+ typename LinearSolver =
+ Eigen::LDLT<Eigen::Matrix<typename Function::Scalar,
+ Function::NUM_PARAMETERS,
+ Function::NUM_PARAMETERS>>>
class TinySolver {
public:
+ // This class needs to have an Eigen aligned operator new as it contains
+ // fixed-size Eigen types.
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+
enum {
NUM_RESIDUALS = Function::NUM_RESIDUALS,
NUM_PARAMETERS = Function::NUM_PARAMETERS
@@ -164,7 +168,7 @@
Status status = HIT_MAX_ITERATIONS;
};
- bool Update(const Function& function, const Parameters &x) {
+ bool Update(const Function& function, const Parameters& x) {
if (!function(x.data(), error_.data(), jacobian_.data())) {
return false;
}
diff --git a/include/ceres/tiny_solver_autodiff_function.h b/include/ceres/tiny_solver_autodiff_function.h
index cc93d73..b782f54 100644
--- a/include/ceres/tiny_solver_autodiff_function.h
+++ b/include/ceres/tiny_solver_autodiff_function.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2017 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -37,8 +37,8 @@
#include <memory>
#include <type_traits>
-#include "Eigen/Core"
+#include "Eigen/Core"
#include "ceres/jet.h"
#include "ceres/types.h" // For kImpossibleValue.
@@ -103,12 +103,16 @@
// solver.Solve(f, &x);
//
// WARNING: The cost function adapter is not thread safe.
-template<typename CostFunctor,
- int kNumResiduals,
- int kNumParameters,
- typename T = double>
+template <typename CostFunctor,
+ int kNumResiduals,
+ int kNumParameters,
+ typename T = double>
class TinySolverAutoDiffFunction {
public:
+ // This class needs to have an Eigen aligned operator new as it contains
+ // as a member a Jet type, which itself has a fixed-size Eigen type as member.
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+
TinySolverAutoDiffFunction(const CostFunctor& cost_functor)
: cost_functor_(cost_functor) {
Initialize<kNumResiduals>(cost_functor);
@@ -122,9 +126,7 @@
// This is similar to AutoDifferentiate(), but since there is only one
// parameter block it is easier to inline to avoid overhead.
- bool operator()(const T* parameters,
- T* residuals,
- T* jacobian) const {
+ bool operator()(const T* parameters, T* residuals, T* jacobian) const {
if (jacobian == NULL) {
// No jacobian requested, so just directly call the cost function with
// doubles, skipping jets and derivatives.
@@ -150,9 +152,7 @@
// Copy the jacobian out of the derivative part of the residual jets.
Eigen::Map<Eigen::Matrix<T, kNumResiduals, kNumParameters>> jacobian_matrix(
- jacobian,
- num_residuals_,
- kNumParameters);
+ jacobian, num_residuals_, kNumParameters);
for (int r = 0; r < num_residuals_; ++r) {
residuals[r] = jet_residuals_[r].a;
// Note that while this looks like a fast vectorized write, in practice it
@@ -186,7 +186,7 @@
// The number of residuals is dynamically sized and the number of
// parameters is statically sized.
- template<int R>
+ template <int R>
typename std::enable_if<(R == Eigen::Dynamic), void>::type Initialize(
const CostFunctor& function) {
jet_residuals_.resize(function.NumResiduals());
@@ -194,7 +194,7 @@
}
// The number of parameters and residuals are statically sized.
- template<int R>
+ template <int R>
typename std::enable_if<(R != Eigen::Dynamic), void>::type Initialize(
const CostFunctor& /* function */) {
num_residuals_ = kNumResiduals;
diff --git a/include/ceres/tiny_solver_cost_function_adapter.h b/include/ceres/tiny_solver_cost_function_adapter.h
index d44bdeb..18ccb39 100644
--- a/include/ceres/tiny_solver_cost_function_adapter.h
+++ b/include/ceres/tiny_solver_cost_function_adapter.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2017 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -28,7 +28,6 @@
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
-
#ifndef CERES_PUBLIC_TINY_SOLVER_COST_FUNCTION_ADAPTER_H_
#define CERES_PUBLIC_TINY_SOLVER_COST_FUNCTION_ADAPTER_H_
@@ -72,7 +71,8 @@
//
// TinySolverCostFunctionAdapter cost_function_adapter(*cost_function);
//
-template <int kNumResiduals = Eigen::Dynamic, int kNumParameters = Eigen::Dynamic>
+template <int kNumResiduals = Eigen::Dynamic,
+ int kNumParameters = Eigen::Dynamic>
class TinySolverCostFunctionAdapter {
public:
typedef double Scalar;
@@ -81,6 +81,10 @@
NUM_RESIDUALS = kNumResiduals
};
+ // This struct needs to have an Eigen aligned operator new as it contains
+ // fixed-size Eigen types.
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+
TinySolverCostFunctionAdapter(const CostFunction& cost_function)
: cost_function_(cost_function) {
CHECK_EQ(cost_function_.parameter_block_sizes().size(), 1)
@@ -127,6 +131,7 @@
return cost_function_.parameter_block_sizes()[0];
}
+ private:
const CostFunction& cost_function_;
mutable Eigen::Matrix<double, NUM_RESIDUALS, NUM_PARAMETERS, Eigen::RowMajor>
row_major_jacobian_;
diff --git a/include/ceres/types.h b/include/ceres/types.h
index 8e5da6a..5ee6fdc 100644
--- a/include/ceres/types.h
+++ b/include/ceres/types.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -39,8 +39,8 @@
#include <string>
-#include "ceres/internal/port.h"
#include "ceres/internal/disable_warnings.h"
+#include "ceres/internal/port.h"
namespace ceres {
@@ -50,7 +50,7 @@
// delete on it upon completion.
enum Ownership {
DO_NOT_TAKE_OWNERSHIP,
- TAKE_OWNERSHIP
+ TAKE_OWNERSHIP,
};
// TODO(keir): Considerably expand the explanations of each solver type.
@@ -112,10 +112,29 @@
// the scene to determine the sparsity structure of the
// preconditioner. This is done using a clustering algorithm. The
// available visibility clustering algorithms are described below.
- //
- // Note: Requires SuiteSparse.
CLUSTER_JACOBI,
- CLUSTER_TRIDIAGONAL
+ CLUSTER_TRIDIAGONAL,
+
+ // Subset preconditioner is a general purpose preconditioner
+ // linear least squares problems. Given a set of residual blocks,
+ // it uses the corresponding subset of the rows of the Jacobian to
+ // construct a preconditioner.
+ //
+ // Suppose the Jacobian J has been horizontally partitioned as
+ //
+ // J = [P]
+ // [Q]
+ //
+ // Where, Q is the set of rows corresponding to the residual
+ // blocks in residual_blocks_for_subset_preconditioner.
+ //
+ // The preconditioner is the inverse of the matrix Q'Q.
+ //
+ // Obviously, the efficacy of the preconditioner depends on how
+ // well the matrix Q approximates J'J, or how well the chosen
+ // residual blocks approximate the non-linear least squares
+ // problem.
+ SUBSET,
};
enum VisibilityClusteringType {
@@ -166,19 +185,19 @@
enum DenseLinearAlgebraLibraryType {
EIGEN,
- LAPACK
+ LAPACK,
};
// Logging options
// The options get progressively noisier.
enum LoggingType {
SILENT,
- PER_MINIMIZER_ITERATION
+ PER_MINIMIZER_ITERATION,
};
enum MinimizerType {
LINE_SEARCH,
- TRUST_REGION
+ TRUST_REGION,
};
enum LineSearchDirectionType {
@@ -221,26 +240,26 @@
// For more details on BFGS see:
//
// Broyden, C.G., "The Convergence of a Class of Double-rank Minimization
- // Algorithms,"; J. Inst. Maths. Applics., Vol. 6, pp 76–90, 1970.
+ // Algorithms,"; J. Inst. Maths. Applics., Vol. 6, pp 76-90, 1970.
//
// Fletcher, R., "A New Approach to Variable Metric Algorithms,"
- // Computer Journal, Vol. 13, pp 317–322, 1970.
+ // Computer Journal, Vol. 13, pp 317-322, 1970.
//
// Goldfarb, D., "A Family of Variable Metric Updates Derived by Variational
- // Means," Mathematics of Computing, Vol. 24, pp 23–26, 1970.
+ // Means," Mathematics of Computing, Vol. 24, pp 23-26, 1970.
//
// Shanno, D.F., "Conditioning of Quasi-Newton Methods for Function
- // Minimization," Mathematics of Computing, Vol. 24, pp 647–656, 1970.
+ // Minimization," Mathematics of Computing, Vol. 24, pp 647-656, 1970.
//
// For more details on L-BFGS see:
//
// Nocedal, J. (1980). "Updating Quasi-Newton Matrices with Limited
- // Storage". Mathematics of Computation 35 (151): 773–782.
+ // Storage". Mathematics of Computation 35 (151): 773-782.
//
// Byrd, R. H.; Nocedal, J.; Schnabel, R. B. (1994).
// "Representations of Quasi-Newton Matrices and their use in
// Limited Memory Methods". Mathematical Programming 63 (4):
- // 129–156.
+ // 129-156.
//
// A general reference for both methods:
//
@@ -393,7 +412,7 @@
// specified for the number of residuals. If specified, then the
// number of residuas for that cost function can vary at runtime.
enum DimensionType {
- DYNAMIC = -1
+ DYNAMIC = -1,
};
// The differentiation method used to compute numerical derivatives in
@@ -414,7 +433,7 @@
enum LineSearchInterpolationType {
BISECTION,
QUADRATIC,
- CUBIC
+ CUBIC,
};
enum CovarianceAlgorithmType {
@@ -429,8 +448,7 @@
// did not write to that memory location.
const double kImpossibleValue = 1e302;
-CERES_EXPORT const char* LinearSolverTypeToString(
- LinearSolverType type);
+CERES_EXPORT const char* LinearSolverTypeToString(LinearSolverType type);
CERES_EXPORT bool StringToLinearSolverType(std::string value,
LinearSolverType* type);
@@ -440,25 +458,23 @@
CERES_EXPORT const char* VisibilityClusteringTypeToString(
VisibilityClusteringType type);
-CERES_EXPORT bool StringToVisibilityClusteringType(std::string value,
- VisibilityClusteringType* type);
+CERES_EXPORT bool StringToVisibilityClusteringType(
+ std::string value, VisibilityClusteringType* type);
CERES_EXPORT const char* SparseLinearAlgebraLibraryTypeToString(
SparseLinearAlgebraLibraryType type);
CERES_EXPORT bool StringToSparseLinearAlgebraLibraryType(
- std::string value,
- SparseLinearAlgebraLibraryType* type);
+ std::string value, SparseLinearAlgebraLibraryType* type);
CERES_EXPORT const char* DenseLinearAlgebraLibraryTypeToString(
DenseLinearAlgebraLibraryType type);
CERES_EXPORT bool StringToDenseLinearAlgebraLibraryType(
- std::string value,
- DenseLinearAlgebraLibraryType* type);
+ std::string value, DenseLinearAlgebraLibraryType* type);
CERES_EXPORT const char* TrustRegionStrategyTypeToString(
TrustRegionStrategyType type);
-CERES_EXPORT bool StringToTrustRegionStrategyType(std::string value,
- TrustRegionStrategyType* type);
+CERES_EXPORT bool StringToTrustRegionStrategyType(
+ std::string value, TrustRegionStrategyType* type);
CERES_EXPORT const char* DoglegTypeToString(DoglegType type);
CERES_EXPORT bool StringToDoglegType(std::string value, DoglegType* type);
@@ -468,35 +484,40 @@
CERES_EXPORT const char* LineSearchDirectionTypeToString(
LineSearchDirectionType type);
-CERES_EXPORT bool StringToLineSearchDirectionType(std::string value,
- LineSearchDirectionType* type);
+CERES_EXPORT bool StringToLineSearchDirectionType(
+ std::string value, LineSearchDirectionType* type);
CERES_EXPORT const char* LineSearchTypeToString(LineSearchType type);
-CERES_EXPORT bool StringToLineSearchType(std::string value, LineSearchType* type);
+CERES_EXPORT bool StringToLineSearchType(std::string value,
+ LineSearchType* type);
CERES_EXPORT const char* NonlinearConjugateGradientTypeToString(
NonlinearConjugateGradientType type);
CERES_EXPORT bool StringToNonlinearConjugateGradientType(
- std::string value,
- NonlinearConjugateGradientType* type);
+ std::string value, NonlinearConjugateGradientType* type);
CERES_EXPORT const char* LineSearchInterpolationTypeToString(
LineSearchInterpolationType type);
CERES_EXPORT bool StringToLineSearchInterpolationType(
- std::string value,
- LineSearchInterpolationType* type);
+ std::string value, LineSearchInterpolationType* type);
CERES_EXPORT const char* CovarianceAlgorithmTypeToString(
CovarianceAlgorithmType type);
CERES_EXPORT bool StringToCovarianceAlgorithmType(
- std::string value,
- CovarianceAlgorithmType* type);
+ std::string value, CovarianceAlgorithmType* type);
CERES_EXPORT const char* NumericDiffMethodTypeToString(
NumericDiffMethodType type);
-CERES_EXPORT bool StringToNumericDiffMethodType(
- std::string value,
- NumericDiffMethodType* type);
+CERES_EXPORT bool StringToNumericDiffMethodType(std::string value,
+ NumericDiffMethodType* type);
+
+CERES_EXPORT const char* LoggingTypeToString(LoggingType type);
+CERES_EXPORT bool StringtoLoggingType(std::string value, LoggingType* type);
+
+CERES_EXPORT const char* DumpFormatTypeToString(DumpFormatType type);
+CERES_EXPORT bool StringtoDumpFormatType(std::string value,
+ DumpFormatType* type);
+CERES_EXPORT bool StringtoDumpFormatType(std::string value, LoggingType* type);
CERES_EXPORT const char* TerminationTypeToString(TerminationType type);
diff --git a/include/ceres/version.h b/include/ceres/version.h
index 20fbe22..a76cc10 100644
--- a/include/ceres/version.h
+++ b/include/ceres/version.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2019 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -41,8 +41,9 @@
#define CERES_TO_STRING(x) CERES_TO_STRING_HELPER(x)
// The Ceres version as a string; for example "1.9.0".
-#define CERES_VERSION_STRING CERES_TO_STRING(CERES_VERSION_MAJOR) "." \
- CERES_TO_STRING(CERES_VERSION_MINOR) "." \
- CERES_TO_STRING(CERES_VERSION_REVISION)
+#define CERES_VERSION_STRING \
+ CERES_TO_STRING(CERES_VERSION_MAJOR) \
+ "." CERES_TO_STRING(CERES_VERSION_MINOR) "." CERES_TO_STRING( \
+ CERES_VERSION_REVISION)
#endif // CERES_PUBLIC_VERSION_H_