blender
/
blender
134
400
Fork 583
Code Issues 6.3k Pull Requests 759 Projects 19 Wiki Activity

Update Ceres to latest upstream version 

As usual brings fixes and speed improvements.
This commit is contained in:
Sergey Sharybin 2014-09-24 17:10:02 +06:00
parent 181e7f98b2
commit 42abfe4853
92 changed files with 5410 additions and 3854 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
extern/libmv/third_party/ceres/CMakeLists.txt vendored
View File

@ -51,6 +51,7 @@ set(SRC
internal/ceres/block_random_access_sparse_matrix.cc
internal/ceres/block_sparse_matrix.cc
internal/ceres/block_structure.cc
internal/ceres/callbacks.cc
internal/ceres/canonical_views_clustering.cc
internal/ceres/c_api.cc
internal/ceres/cgnr_solver.cc
@ -76,6 +77,8 @@ set(SRC
internal/ceres/generated/partitioned_matrix_view_d_d_d.cc
internal/ceres/generated/schur_eliminator_d_d_d.cc
internal/ceres/gradient_checking_cost_function.cc
internal/ceres/gradient_problem.cc
internal/ceres/gradient_problem_solver.cc
internal/ceres/implicit_schur_complement.cc
internal/ceres/incomplete_lq_factorization.cc
internal/ceres/iterative_schur_complement_solver.cc
@ -87,6 +90,7 @@ set(SRC
internal/ceres/line_search.cc
internal/ceres/line_search_direction.cc
internal/ceres/line_search_minimizer.cc
internal/ceres/line_search_preprocessor.cc
internal/ceres/local_parameterization.cc
internal/ceres/loss_function.cc
internal/ceres/low_rank_inverse_hessian.cc
@ -96,9 +100,11 @@ set(SRC
internal/ceres/partitioned_matrix_view.cc
internal/ceres/polynomial.cc
internal/ceres/preconditioner.cc
internal/ceres/preprocessor.cc
internal/ceres/problem.cc
internal/ceres/problem_impl.cc
internal/ceres/program.cc
internal/ceres/reorder_program.cc
internal/ceres/residual_block.cc
internal/ceres/residual_block_utils.cc
internal/ceres/schur_complement_solver.cc
@ -107,7 +113,7 @@ set(SRC
internal/ceres/scratch_evaluate_preparer.cc
internal/ceres/single_linkage_clustering.cc
internal/ceres/solver.cc
internal/ceres/solver_impl.cc
internal/ceres/solver_utils.cc
internal/ceres/sparse_matrix.cc
internal/ceres/sparse_normal_cholesky_solver.cc
internal/ceres/split.cc
@ -115,6 +121,7 @@ set(SRC
internal/ceres/suitesparse.cc
internal/ceres/triplet_sparse_matrix.cc
internal/ceres/trust_region_minimizer.cc
internal/ceres/trust_region_preprocessor.cc
internal/ceres/trust_region_strategy.cc
internal/ceres/types.cc
internal/ceres/visibility_based_preconditioner.cc
@ -134,13 +141,17 @@ set(SRC
include/ceres/dynamic_numeric_diff_cost_function.h
include/ceres/fpclassify.h
include/ceres/gradient_checker.h
include/ceres/gradient_problem.h
include/ceres/gradient_problem_solver.h
include/ceres/internal/autodiff.h
include/ceres/internal/disable_warnings.h
include/ceres/internal/eigen.h
include/ceres/internal/fixed_array.h
include/ceres/internal/macros.h
include/ceres/internal/manual_constructor.h
include/ceres/internal/numeric_diff.h
include/ceres/internal/port.h
include/ceres/internal/reenable_warnings.h
include/ceres/internal/scoped_ptr.h
include/ceres/internal/variadic_evaluate.h
include/ceres/iteration_callback.h
@ -149,13 +160,13 @@ set(SRC
include/ceres/loss_function.h
include/ceres/normal_prior.h
include/ceres/numeric_diff_cost_function.h
include/ceres/numeric_diff_functor.h
include/ceres/ordered_groups.h
include/ceres/problem.h
include/ceres/rotation.h
include/ceres/sized_cost_function.h
include/ceres/solver.h
include/ceres/types.h
include/ceres/version.h
internal/ceres/array_utils.h
internal/ceres/blas.h
internal/ceres/block_evaluate_preparer.h
@ -167,6 +178,7 @@ set(SRC
internal/ceres/block_random_access_sparse_matrix.h
internal/ceres/block_sparse_matrix.h
internal/ceres/block_structure.h
internal/ceres/callbacks.h
internal/ceres/canonical_views_clustering.h
internal/ceres/casts.h
internal/ceres/cgnr_linear_operator.h
@ -193,6 +205,7 @@ set(SRC
internal/ceres/execution_summary.h
internal/ceres/file.h
internal/ceres/gradient_checking_cost_function.h
internal/ceres/gradient_problem_evaluator.h
internal/ceres/graph_algorithms.h
internal/ceres/graph.h
internal/ceres/implicit_schur_complement.h
@ -207,6 +220,7 @@ set(SRC
internal/ceres/line_search_direction.h
internal/ceres/line_search.h
internal/ceres/line_search_minimizer.h
internal/ceres/line_search_preprocessor.h
internal/ceres/low_rank_inverse_hessian.h
internal/ceres/map_util.h
internal/ceres/minimizer.h
@ -217,10 +231,12 @@ set(SRC
internal/ceres/partitioned_matrix_view_impl.h
internal/ceres/polynomial.h
internal/ceres/preconditioner.h
internal/ceres/preprocessor.h
internal/ceres/problem_impl.h
internal/ceres/program_evaluator.h
internal/ceres/program.h
internal/ceres/random.h
internal/ceres/reorder_program.h
internal/ceres/residual_block.h
internal/ceres/residual_block_utils.h
internal/ceres/schur_complement_solver.h
@ -230,7 +246,7 @@ set(SRC
internal/ceres/scratch_evaluate_preparer.h
internal/ceres/single_linkage_clustering.h
internal/ceres/small_blas.h
internal/ceres/solver_impl.h
internal/ceres/solver_utils.h
internal/ceres/sparse_matrix.h
internal/ceres/sparse_normal_cholesky_solver.h
internal/ceres/split.h
@ -239,6 +255,7 @@ set(SRC
internal/ceres/suitesparse.h
internal/ceres/triplet_sparse_matrix.h
internal/ceres/trust_region_minimizer.h
internal/ceres/trust_region_preprocessor.h
internal/ceres/trust_region_strategy.h
internal/ceres/visibility_based_preconditioner.h
internal/ceres/visibility.h

1094
extern/libmv/third_party/ceres/ChangeLog vendored
View File

File diff suppressed because it is too large Load Diff

10
extern/libmv/third_party/ceres/bundle.sh vendored
View File

@ -46,7 +46,7 @@ rm -rf $tmp
sources=`find ./include ./internal -type f -iname '*.cc' -or -iname '*.cpp' -or -iname '*.c' | sed -r 's/^\.\//\t/' | \
grep -v -E 'schur_eliminator_[0-9]_[0-9d]_[0-9d].cc' | \
grep -v -E 'partitioned_matrix_view_[0-9]_[0-9d]_[0-9d].cc' | sort -d`
generated_sources=`find ./include ./internal -type f -iname '*.cc' -or -iname '*.cpp' -or -iname '*.c' | sed -r 's/^\.\//#\t\t/' | \
generated_sources=`find ./include ./internal -type f -iname '*.cc' -or -iname '*.cpp' -or -iname '*.c' | sed -r 's/^\.\//\t\t/' | \
grep -E 'schur_eliminator_[0-9]_[0-9d]_[0-9d].cc|partitioned_matrix_view_[0-9]_[0-9d]_[0-9d].cc' | sort -d`
headers=`find ./include ./internal -type f -iname '*.h' | sed -r 's/^\.\//\t/' | sort -d`
@ -138,11 +138,11 @@ ${sources}
${headers}
)
#if(TRUE)
# list(APPEND SRC
if(TRUE)
list(APPEND SRC
${generated_sources}
# )
#endif()
)
endif()
if(WIN32)
list(APPEND INC

24
extern/libmv/third_party/ceres/files.txt vendored
View File

@ -11,13 +11,17 @@ include/ceres/dynamic_autodiff_cost_function.h
include/ceres/dynamic_numeric_diff_cost_function.h
include/ceres/fpclassify.h
include/ceres/gradient_checker.h
include/ceres/gradient_problem.h
include/ceres/gradient_problem_solver.h
include/ceres/internal/autodiff.h
include/ceres/internal/disable_warnings.h
include/ceres/internal/eigen.h
include/ceres/internal/fixed_array.h
include/ceres/internal/macros.h
include/ceres/internal/manual_constructor.h
include/ceres/internal/numeric_diff.h
include/ceres/internal/port.h
include/ceres/internal/reenable_warnings.h
include/ceres/internal/scoped_ptr.h
include/ceres/internal/variadic_evaluate.h
include/ceres/iteration_callback.h
@ -26,13 +30,13 @@ include/ceres/local_parameterization.h
include/ceres/loss_function.h
include/ceres/normal_prior.h
include/ceres/numeric_diff_cost_function.h
include/ceres/numeric_diff_functor.h
include/ceres/ordered_groups.h
include/ceres/problem.h
include/ceres/rotation.h
include/ceres/sized_cost_function.h
include/ceres/solver.h
include/ceres/types.h
include/ceres/version.h
internal/ceres/array_utils.cc
internal/ceres/array_utils.h
internal/ceres/blas.cc
@ -55,6 +59,8 @@ internal/ceres/block_sparse_matrix.cc
internal/ceres/block_sparse_matrix.h
internal/ceres/block_structure.cc
internal/ceres/block_structure.h
internal/ceres/callbacks.cc
internal/ceres/callbacks.h
internal/ceres/canonical_views_clustering.cc
internal/ceres/canonical_views_clustering.h
internal/ceres/c_api.cc
@ -62,7 +68,6 @@ internal/ceres/casts.h
internal/ceres/cgnr_linear_operator.h
internal/ceres/cgnr_solver.cc
internal/ceres/cgnr_solver.h
internal/ceres/CMakeLists.txt
internal/ceres/collections_port.h
internal/ceres/compressed_col_sparse_matrix_utils.cc
internal/ceres/compressed_col_sparse_matrix_utils.h
@ -145,6 +150,9 @@ internal/ceres/generate_eliminator_specialization.py
internal/ceres/generate_partitioned_matrix_view_specializations.py
internal/ceres/gradient_checking_cost_function.cc
internal/ceres/gradient_checking_cost_function.h
internal/ceres/gradient_problem.cc
internal/ceres/gradient_problem_evaluator.h
internal/ceres/gradient_problem_solver.cc
internal/ceres/graph_algorithms.h
internal/ceres/graph.h
internal/ceres/implicit_schur_complement.cc
@ -170,6 +178,8 @@ internal/ceres/line_search_direction.h
internal/ceres/line_search.h
internal/ceres/line_search_minimizer.cc
internal/ceres/line_search_minimizer.h
internal/ceres/line_search_preprocessor.cc
internal/ceres/line_search_preprocessor.h
internal/ceres/local_parameterization.cc
internal/ceres/loss_function.cc
internal/ceres/low_rank_inverse_hessian.cc
@ -189,6 +199,8 @@ internal/ceres/polynomial.cc
internal/ceres/polynomial.h
internal/ceres/preconditioner.cc
internal/ceres/preconditioner.h
internal/ceres/preprocessor.cc
internal/ceres/preprocessor.h
internal/ceres/problem.cc
internal/ceres/problem_impl.cc
internal/ceres/problem_impl.h
@ -196,6 +208,8 @@ internal/ceres/program.cc
internal/ceres/program_evaluator.h
internal/ceres/program.h
internal/ceres/random.h
internal/ceres/reorder_program.cc
internal/ceres/reorder_program.h
internal/ceres/residual_block.cc
internal/ceres/residual_block.h
internal/ceres/residual_block_utils.cc
@ -213,8 +227,8 @@ internal/ceres/single_linkage_clustering.cc
internal/ceres/single_linkage_clustering.h
internal/ceres/small_blas.h
internal/ceres/solver.cc
internal/ceres/solver_impl.cc
internal/ceres/solver_impl.h
internal/ceres/solver_utils.cc
internal/ceres/solver_utils.h
internal/ceres/sparse_matrix.cc
internal/ceres/sparse_matrix.h
internal/ceres/sparse_normal_cholesky_solver.cc
@ -230,6 +244,8 @@ internal/ceres/triplet_sparse_matrix.cc
internal/ceres/triplet_sparse_matrix.h
internal/ceres/trust_region_minimizer.cc
internal/ceres/trust_region_minimizer.h
internal/ceres/trust_region_preprocessor.cc
internal/ceres/trust_region_preprocessor.h
internal/ceres/trust_region_strategy.cc
internal/ceres/trust_region_strategy.h
internal/ceres/types.cc

3
extern/libmv/third_party/ceres/include/ceres/c_api.h vendored
View File

@ -39,6 +39,7 @@
#define CERES_PUBLIC_C_API_H_
#include "ceres/internal/port.h"
#include "ceres/internal/disable_warnings.h"
#ifdef __cplusplus
extern "C" {
@ -140,4 +141,6 @@ CERES_EXPORT void ceres_solve(ceres_problem_t* problem);
}
#endif
#include "ceres/internal/reenable_warnings.h"
#endif /* CERES_PUBLIC_C_API_H_ */

9
extern/libmv/third_party/ceres/include/ceres/ceres.h vendored
View File

@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
@ -34,9 +34,6 @@
#ifndef CERES_PUBLIC_CERES_H_
#define CERES_PUBLIC_CERES_H_
#define CERES_VERSION 1.9.0
#define CERES_ABI_VERSION 1.9.0
#include "ceres/autodiff_cost_function.h"
#include "ceres/autodiff_local_parameterization.h"
#include "ceres/cost_function.h"
@ -45,16 +42,18 @@
#include "ceres/crs_matrix.h"
#include "ceres/dynamic_autodiff_cost_function.h"
#include "ceres/dynamic_numeric_diff_cost_function.h"
#include "ceres/gradient_problem.h"
#include "ceres/gradient_problem_solver.h"
#include "ceres/iteration_callback.h"
#include "ceres/jet.h"
#include "ceres/local_parameterization.h"
#include "ceres/loss_function.h"
#include "ceres/numeric_diff_cost_function.h"
#include "ceres/numeric_diff_functor.h"
#include "ceres/ordered_groups.h"
#include "ceres/problem.h"
#include "ceres/sized_cost_function.h"
#include "ceres/solver.h"
#include "ceres/types.h"
#include "ceres/version.h"
#endif // CERES_PUBLIC_CERES_H_

2
extern/libmv/third_party/ceres/include/ceres/conditioned_cost_function.h vendored
View File

@ -39,6 +39,7 @@
#include "ceres/cost_function.h"
#include "ceres/internal/scoped_ptr.h"
#include "ceres/types.h"
#include "ceres/internal/disable_warnings.h"
namespace ceres {
@ -93,5 +94,6 @@ class CERES_EXPORT ConditionedCostFunction : public CostFunction {
} // namespace ceres
#include "ceres/internal/reenable_warnings.h"
#endif // CERES_PUBLIC_CONDITIONED_COST_FUNCTION_H_

6
extern/libmv/third_party/ceres/include/ceres/cost_function.h vendored
View File

@ -48,6 +48,7 @@
#include "ceres/internal/macros.h"
#include "ceres/internal/port.h"
#include "ceres/types.h"
#include "ceres/internal/disable_warnings.h"
namespace ceres {
@ -105,8 +106,7 @@ class CERES_EXPORT CostFunction {
// the constraints, then returning false whenever the constraints
// are not satisfied will prevent the solver from moving into the
// infeasible region. This is not a very sophisticated mechanism for
// enforcing constraints, but is often good enough for things like
// non-negativity constraints.
// enforcing constraints, but is often good enough.
//
// Note that it is important that the initial values of the
// parameter block must be feasible, otherwise the solver will
@ -142,4 +142,6 @@ class CERES_EXPORT CostFunction {
} // namespace ceres
#include "ceres/internal/reenable_warnings.h"
#endif // CERES_PUBLIC_COST_FUNCTION_H_

4
extern/libmv/third_party/ceres/include/ceres/cost_function_to_functor.h vendored
View File

@ -107,9 +107,7 @@ class CostFunctionToFunctor {
explicit CostFunctionToFunctor(CostFunction* cost_function)
: cost_function_(cost_function) {
CHECK_NOTNULL(cost_function);
CHECK_GE(kNumResiduals, 0);
CHECK_EQ(cost_function->num_residuals(), kNumResiduals);
CHECK(kNumResiduals > 0 || kNumResiduals == DYNAMIC);
// This block breaks the 80 column rule to keep it somewhat readable.
CHECK((!N1 && !N2 && !N3 && !N4 && !N5 && !N6 && !N7 && !N8 && !N9) ||

78
extern/libmv/third_party/ceres/include/ceres/covariance.h vendored
View File

@ -36,6 +36,7 @@
#include "ceres/internal/port.h"
#include "ceres/internal/scoped_ptr.h"
#include "ceres/types.h"
#include "ceres/internal/disable_warnings.h"
namespace ceres {
@ -201,9 +202,9 @@ class CERES_EXPORT Covariance {
struct CERES_EXPORT Options {
Options()
#ifndef CERES_NO_SUITESPARSE
: algorithm_type(SPARSE_QR),
: algorithm_type(SUITE_SPARSE_QR),
#else
: algorithm_type(DENSE_SVD),
: algorithm_type(EIGEN_SPARSE_QR),
#endif
min_reciprocal_condition_number(1e-14),
null_space_rank(0),
@ -228,47 +229,22 @@ class CERES_EXPORT Covariance {
// for small to moderate sized problems. It can handle
// full-rank as well as rank deficient Jacobians.
//
// 2. SPARSE_CHOLESKY uses the CHOLMOD sparse Cholesky
// factorization library to compute the decomposition :
//
// R'R = J'J
//
// and then
//
// [J'J]^-1 = [R'R]^-1
//
// It a fast algorithm for sparse matrices that should be used
// when the Jacobian matrix J is well conditioned. For
// ill-conditioned matrices, this algorithm can fail
// unpredictabily. This is because Cholesky factorization is
// not a rank-revealing factorization, i.e., it cannot reliably
// detect when the matrix being factorized is not of full
// rank. SuiteSparse/CHOLMOD supplies a heuristic for checking
// if the matrix is rank deficient (cholmod_rcond), but it is
// only a heuristic and can have both false positive and false
// negatives.
//
// Recent versions of SuiteSparse (>= 4.2.0) provide a much
// more efficient method for solving for rows of the covariance
// matrix. Therefore, if you are doing SPARSE_CHOLESKY, we
// strongly recommend using a recent version of SuiteSparse.
//
// 3. SPARSE_QR uses the SuiteSparseQR sparse QR factorization
// library to compute the decomposition
// 2. EIGEN_SPARSE_QR uses the sparse QR factorization algorithm
// in Eigen to compute the decomposition
//
// Q * R = J
//
// [J'J]^-1 = [R*R']^-1
//
// It is a moderately fast algorithm for sparse matrices, which
// at the price of more time and memory than the
// SPARSE_CHOLESKY algorithm is numerically better behaved and
// is rank revealing, i.e., it can reliably detect when the
// Jacobian matrix is rank deficient.
// It is a moderately fast algorithm for sparse matrices.
//
// Neither SPARSE_CHOLESKY or SPARSE_QR are capable of computing
// the covariance if the Jacobian is rank deficient.
// 3. SUITE_SPARSE_QR uses the SuiteSparseQR sparse QR
// factorization algorithm. It uses dense linear algebra and is
// multi threaded, so for large sparse sparse matrices it is
// significantly faster than EIGEN_SPARSE_QR.
//
// Neither EIGEN_SPARSE_QR not SUITE_SPARSE_QR are capable of
// computing the covariance if the Jacobian is rank deficient.
CovarianceAlgorithmType algorithm_type;
// If the Jacobian matrix is near singular, then inverting J'J
@ -294,29 +270,13 @@ class CERES_EXPORT Covariance {
// where min_sigma and max_sigma are the minimum and maxiumum
// singular values of J respectively.
//
// 2. SPARSE_CHOLESKY
//
// cholmod_rcond < min_reciprocal_conditioner_number
//
// Here cholmod_rcond is a crude estimate of the reciprocal
// condition number of J'J by using the maximum and minimum
// diagonal entries of the Cholesky factor R. There are no
// theoretical guarantees associated with this test. It can
// give false positives and negatives. Use at your own
// risk. The default value of min_reciprocal_condition_number
// has been set to a conservative value, and sometimes the
// Covariance::Compute may return false even if it is possible
// to estimate the covariance reliably. In such cases, the user
// should exercise their judgement before lowering the value of
// min_reciprocal_condition_number.
//
// 3. SPARSE_QR
// 2. SUITE_SPARSE_QR and EIGEN_SPARSE_QR
//
// rank(J) < num_col(J)
//
// Here rank(J) is the estimate of the rank of J returned by the
// SuiteSparseQR algorithm. It is a fairly reliable indication
// of rank deficiency.
// sparse QR factorization algorithm. It is a fairly reliable
// indication of rank deficiency.
//
double min_reciprocal_condition_number;
@ -351,8 +311,8 @@ class CERES_EXPORT Covariance {
//
// lambda_i / lambda_max < min_reciprocal_condition_number.
//
// This option has no effect on the SPARSE_CHOLESKY or SPARSE_QR
// algorithms.
// This option has no effect on the SUITE_SPARSE_QR and
// EIGEN_SPARSE_QR algorithms.
int null_space_rank;
int num_threads;
@ -419,4 +379,6 @@ class CERES_EXPORT Covariance {
} // namespace ceres
#include "ceres/internal/reenable_warnings.h"
#endif // CERES_PUBLIC_COVARIANCE_H_

3
extern/libmv/third_party/ceres/include/ceres/crs_matrix.h vendored
View File

@ -33,6 +33,7 @@
#include <vector>
#include "ceres/internal/port.h"
#include "ceres/internal/disable_warnings.h"
namespace ceres {
@ -80,4 +81,6 @@ struct CERES_EXPORT CRSMatrix {
} // namespace ceres
#include "ceres/internal/reenable_warnings.h"
#endif // CERES_PUBLIC_CRS_MATRIX_H_

127
extern/libmv/third_party/ceres/include/ceres/gradient_problem.h vendored Normal file
View File

@ -0,0 +1,127 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#ifndef CERES_PUBLIC_GRADIENT_PROBLEM_H_
#define CERES_PUBLIC_GRADIENT_PROBLEM_H_
#include "ceres/internal/macros.h"
#include "ceres/internal/port.h"
#include "ceres/internal/scoped_ptr.h"
#include "ceres/local_parameterization.h"
namespace ceres {
class FirstOrderFunction;
// Instances of GradientProblem represent general non-linear
// optimization problems that must be solved using just the value of
// the objective function and its gradient. Unlike the Problem class,
// which can only be used to model non-linear least squares problems,
// instances of GradientProblem not restricted in the form of the
// objective function.
//
// Structurally GradientProblem is a composition of a
// FirstOrderFunction and optionally a LocalParameterization.
//
// The FirstOrderFunction is responsible for evaluating the cost and
// gradient of the objective function.
//
// The LocalParameterization is responsible for going back and forth
// between the ambient space and the local tangent space. (See
// local_parameterization.h for more details). When a
// LocalParameterization is not provided, then the tangent space is
// assumed to coincide with the ambient Euclidean space that the
// gradient vector lives in.
//
// Example usage:
//
// The following demonstrate the problem construction for Rosenbrock's function
//
// f(x,y) = (1-x)^2 + 100(y - x^2)^2;
//
// class Rosenbrock : public ceres::FirstOrderFunction {
// public:
// virtual ~Rosenbrock() {}
//
// virtual bool Evaluate(const double* parameters,
// double* cost,
// double* gradient) const {
// const double x = parameters[0];
// const double y = parameters[1];
//
// cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x);
// if (gradient != NULL) {
// gradient[0] = -2.0 * (1.0 - x) - 200.0 * (y - x * x) * 2.0 * x;
// gradient[1] = 200.0 * (y - x * x);
// }
// return true;
// };
//
// virtual int NumParameters() const { return 2; };
// };
//
// ceres::GradientProblem problem(new Rosenbrock());
class CERES_EXPORT GradientProblem {
public:
// Takes ownership of the function.
explicit GradientProblem(FirstOrderFunction* function);
// Takes ownership of the function and the parameterization.
GradientProblem(FirstOrderFunction* function,
LocalParameterization* parameterization);
int NumParameters() const;
int NumLocalParameters() const;
// This call is not thread safe.
bool Evaluate(const double* parameters, double* cost, double* gradient) const;
bool Plus(const double* x, const double* delta, double* x_plus_delta) const;
private:
internal::scoped_ptr<FirstOrderFunction> function_;
internal::scoped_ptr<LocalParameterization> parameterization_;
internal::scoped_array<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_

365
extern/libmv/third_party/ceres/include/ceres/gradient_problem_solver.h vendored Normal file
View File

@ -0,0 +1,365 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#ifndef CERES_PUBLIC_GRADIENT_PROBLEM_SOLVER_H_
#define CERES_PUBLIC_GRADIENT_PROBLEM_SOLVER_H_
#include <cmath>
#include <string>
#include <vector>
#include "ceres/internal/macros.h"
#include "ceres/internal/port.h"
#include "ceres/iteration_callback.h"
#include "ceres/types.h"
#include "ceres/internal/disable_warnings.h"
namespace ceres {
class GradientProblem;
class CERES_EXPORT GradientProblemSolver {
public:
virtual ~GradientProblemSolver();
// The options structure contains, not surprisingly, options that control how
// the solver operates. The defaults should be suitable for a wide range of
// problems; however, better performance is often obtainable with tweaking.
//
// The constants are defined inside types.h
struct CERES_EXPORT Options {
// Default constructor that sets up a generic sparse problem.
Options() {
line_search_direction_type = LBFGS;
line_search_type = WOLFE;
nonlinear_conjugate_gradient_type = FLETCHER_REEVES;
max_lbfgs_rank = 20;
use_approximate_eigenvalue_bfgs_scaling = false;
line_search_interpolation_type = CUBIC;
min_line_search_step_size = 1e-9;
line_search_sufficient_function_decrease = 1e-4;
max_line_search_step_contraction = 1e-3;
min_line_search_step_contraction = 0.6;
max_num_line_search_step_size_iterations = 20;
max_num_line_search_direction_restarts = 5;
line_search_sufficient_curvature_decrease = 0.9;
max_line_search_step_expansion = 10.0;
max_num_iterations = 50;
max_solver_time_in_seconds = 1e9;
num_threads = 1;
function_tolerance = 1e-6;
gradient_tolerance = 1e-10;
logging_type = PER_MINIMIZER_ITERATION;
minimizer_progress_to_stdout = false;
}
// Returns true if the options struct has a valid
// configuration. Returns false otherwise, and fills in *error
// with a message describing the problem.
bool IsValid(string* error) const;
// Minimizer options ----------------------------------------
LineSearchDirectionType line_search_direction_type;
LineSearchType line_search_type;
NonlinearConjugateGradientType nonlinear_conjugate_gradient_type;
// The LBFGS hessian approximation is a low rank approximation to
// the inverse of the Hessian matrix. The rank of the
// approximation determines (linearly) the space and time
// complexity of using the approximation. Higher the rank, the
// better is the quality of the approximation. The increase in
// quality is however is bounded for a number of reasons.
//
// 1. The method only uses secant information and not actual
// derivatives.
//
// 2. The Hessian approximation is constrained to be positive
// definite.
//
// So increasing this rank to a large number will cost time and
// space complexity without the corresponding increase in solution
// quality. There are no hard and fast rules for choosing the
// maximum rank. The best choice usually requires some problem
// specific experimentation.
//
// For more theoretical and implementation details of the LBFGS
// method, please see:
//
// Nocedal, J. (1980). "Updating Quasi-Newton Matrices with
// Limited Storage". Mathematics of Computation 35 (151): 773–782.
int max_lbfgs_rank;
// As part of the (L)BFGS update step (BFGS) / right-multiply step (L-BFGS),
// the initial inverse Hessian approximation is taken to be the Identity.
// However, Oren showed that using instead I * \gamma, where \gamma is
// chosen to approximate an eigenvalue of the true inverse Hessian can
// result in improved convergence in a wide variety of cases. Setting
// use_approximate_eigenvalue_bfgs_scaling to true enables this scaling.
//
// It is important to note that approximate eigenvalue scaling does not
// always improve convergence, and that it can in fact significantly degrade
// performance for certain classes of problem, which is why it is disabled
// by default. In particular it can degrade performance when the
// sensitivity of the problem to different parameters varies significantly,
// as in this case a single scalar factor fails to capture this variation
// and detrimentally downscales parts of the jacobian approximation which
// correspond to low-sensitivity parameters. It can also reduce the
// robustness of the solution to errors in the jacobians.
//
// Oren S.S., Self-scaling variable metric (SSVM) algorithms
// Part II: Implementation and experiments, Management Science,
// 20(5), 863-874, 1974.
bool use_approximate_eigenvalue_bfgs_scaling;
// Degree of the polynomial used to approximate the objective
// function. Valid values are BISECTION, QUADRATIC and CUBIC.
//
// BISECTION corresponds to pure backtracking search with no
// interpolation.
LineSearchInterpolationType line_search_interpolation_type;
// If during the line search, the step_size falls below this
// value, it is truncated to zero.
double min_line_search_step_size;
// Line search parameters.
// Solving the line search problem exactly is computationally
// prohibitive. Fortunately, line search based optimization
// algorithms can still guarantee convergence if instead of an
// exact solution, the line search algorithm returns a solution
// which decreases the value of the objective function
// sufficiently. More precisely, we are looking for a step_size
// s.t.
//
// f(step_size) <= f(0) + sufficient_decrease * f'(0) * step_size
//
double line_search_sufficient_function_decrease;
// In each iteration of the line search,
//
// new_step_size >= max_line_search_step_contraction * step_size
//
// Note that by definition, for contraction:
//
// 0 < max_step_contraction < min_step_contraction < 1
//
double max_line_search_step_contraction;
// In each iteration of the line search,
//
// new_step_size <= min_line_search_step_contraction * step_size
//
// Note that by definition, for contraction:
//
// 0 < max_step_contraction < min_step_contraction < 1
//
double min_line_search_step_contraction;
// 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.
int max_num_line_search_step_size_iterations;
// Maximum number of restarts of the line search direction algorithm before
// terminating the optimization. Restarts of the line search direction
// algorithm occur when the current algorithm fails to produce a new descent
// direction. This typically indicates a numerical failure, or a breakdown
// in the validity of the approximations used.
int max_num_line_search_direction_restarts;
// The strong Wolfe conditions consist of the Armijo sufficient
// decrease condition, and an additional requirement that the
// step-size be chosen s.t. the _magnitude_ ('strong' Wolfe
// conditions) of the gradient along the search direction
// decreases sufficiently. Precisely, this second condition
// is that we seek a step_size s.t.
//
// |f'(step_size)| <= sufficient_curvature_decrease * |f'(0)|
//
// Where f() is the line search objective and f'() is the derivative
// of f w.r.t step_size (d f / d step_size).
double line_search_sufficient_curvature_decrease;
// During the bracketing phase of the Wolfe search, the step size is
// increased until either a point satisfying the Wolfe conditions is
// found, or an upper bound for a bracket containing a point satisfying
// the conditions is found. Precisely, at each iteration of the
// expansion:
//
// new_step_size <= max_step_expansion * step_size.
//
// By definition for expansion, max_step_expansion > 1.0.
double max_line_search_step_expansion;
// Maximum number of iterations for the minimizer to run for.
int max_num_iterations;
// Maximum time for which the minimizer should run for.
double max_solver_time_in_seconds;
// Number of threads used by Ceres for evaluating the cost and
// jacobians.
int num_threads;
// Minimizer terminates when
//
// (new_cost - old_cost) < function_tolerance * old_cost;
//
double function_tolerance;
// Minimizer terminates when
//
// max_i |x - Project(Plus(x, -g(x))| < gradient_tolerance
//
// This value should typically be 1e-4 * function_tolerance.
double gradient_tolerance;
// Logging options ---------------------------------------------------------
LoggingType logging_type;
// By default the Minimizer progress is logged to VLOG(1), which
// is sent to STDERR depending on the vlog level. If this flag is
// set to true, and logging_type is not SILENT, the logging output
// is sent to STDOUT.
bool minimizer_progress_to_stdout;
// 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.
bool update_state_every_iteration;
// Callbacks that are executed at the end of each iteration of the
// Minimizer. An iteration may terminate midway, either due to
// numerical failures or because one of the convergence tests has
// been satisfied. In this case none of the callbacks are
// executed.
// Callbacks are executed in the order that they are specified in
// this vector. By default, parameter blocks are updated only at
// the end of the optimization, i.e when the Minimizer
// terminates. This behaviour is controlled by
// update_state_every_variable. If the user wishes to have access
// to the update parameter blocks when his/her callbacks are
// executed, then set update_state_every_iteration to true.
//
// The solver does NOT take ownership of these pointers.
vector<IterationCallback*> callbacks;
};
struct CERES_EXPORT Summary {
Summary();
// A brief one line description of the state of the solver after
// termination.
string BriefReport() const;
// A full multiline description of the state of the solver after
// termination.
string FullReport() const;
bool IsSolutionUsable() const;
// Minimizer summary -------------------------------------------------
TerminationType termination_type;
// Reason why the solver terminated.
string message;
// Cost of the problem (value of the objective function) before
// the optimization.
double initial_cost;
// Cost of the problem (value of the objective function) after the
// optimization.
double final_cost;
// IterationSummary for each minimizer iteration in order.
vector<IterationSummary> iterations;
// Sum total of all time spent inside Ceres when Solve is called.
double total_time_in_seconds;
// Time (in seconds) spent evaluating the residual vector.
double cost_evaluation_time_in_seconds;
// Time (in seconds) spent evaluating the jacobian matrix.
double gradient_evaluation_time_in_seconds;
// Number of parameters in the probem.
int num_parameters;
// Dimension of the tangent space of the problem.
int num_local_parameters;
// Type of line search direction used.
LineSearchDirectionType line_search_direction_type;
// Type of the line search algorithm used.
LineSearchType line_search_type;
// When performing line search, the degree of the polynomial used
// to approximate the objective function.
LineSearchInterpolationType line_search_interpolation_type;
// 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;
// If the type of the line search direction is LBFGS, then this
// indicates the rank of the Hessian approximation.
int max_lbfgs_rank;
};
// Once a least squares problem has been built, this function takes
// the problem and optimizes it based on the values of the options
// parameters. Upon return, a detailed summary of the work performed
// by the preprocessor, the non-linear minmizer and the linear
// solver are reported in the summary object.
virtual void Solve(const GradientProblemSolver::Options& options,
const GradientProblem& problem,
double* parameters,
GradientProblemSolver::Summary* summary);
};
// Helper function which avoids going through the interface.
CERES_EXPORT void Solve(const GradientProblemSolver::Options& options,
const GradientProblem& problem,
double* parameters,
GradientProblemSolver::Summary* summary);
} // namespace ceres
#include "ceres/internal/reenable_warnings.h"
#endif // CERES_PUBLIC_GRADIENT_PROBLEM_SOLVER_H_

44
extern/libmv/third_party/ceres/include/ceres/internal/disable_warnings.h vendored Normal file
View File

@ -0,0 +1,44 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 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.
//
// This file has the sole purpose to silence warnings when including Ceres.
// This is not your usual header guard. The macro CERES_WARNINGS_DISABLED
// shows up again in reenable_warnings.h.
#ifndef CERES_WARNINGS_DISABLED
#define CERES_WARNINGS_DISABLED
#ifdef _MSC_VER
#pragma warning( push )
// Disable the warning C4251 which is trigerred 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 )
#endif
#endif // CERES_WARNINGS_DISABLED

9
extern/libmv/third_party/ceres/include/ceres/internal/numeric_diff.h vendored
View File

@ -103,14 +103,17 @@ struct NumericDiff {
typedef Matrix<double, kNumResiduals, 1> ResidualVector;
typedef Matrix<double, kParameterBlockSize, 1> ParameterVector;
// The convoluted reasoning for choosing the Row/Column major
// ordering of the matrix is an artifact of the restrictions in
// Eigen that prevent it from creating RowMajor matrices with a
// single column. In these cases, we ask for a ColMajor matrix.
typedef Matrix<double,
kNumResiduals,
kParameterBlockSize,
(kParameterBlockSize == 1 &&
kNumResiduals > 1) ? ColMajor : RowMajor>
(kParameterBlockSize == 1) ? ColMajor : RowMajor>
JacobianMatrix;
Map<JacobianMatrix> parameter_jacobian(jacobian,
NUM_RESIDUALS,
kParameterBlockSize);

38
extern/libmv/third_party/ceres/include/ceres/internal/reenable_warnings.h vendored Normal file
View File

@ -0,0 +1,38 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 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.
//
// This is not your usual header guard. See disable_warnings.h
#ifdef CERES_WARNINGS_DISABLED
#undef CERES_WARNINGS_DISABLED
#ifdef _MSC_VER
#pragma warning( pop )
#endif
#endif // CERES_WARNINGS_DISABLED

3
extern/libmv/third_party/ceres/include/ceres/iteration_callback.h vendored
View File

@ -36,6 +36,7 @@
#define CERES_PUBLIC_ITERATION_CALLBACK_H_
#include "ceres/types.h"
#include "ceres/internal/disable_warnings.h"
namespace ceres {
@ -219,4 +220,6 @@ class CERES_EXPORT IterationCallback {
} // namespace ceres
#include "ceres/internal/reenable_warnings.h"
#endif // CERES_PUBLIC_ITERATION_CALLBACK_H_

11
extern/libmv/third_party/ceres/include/ceres/jet.h vendored
View File

@ -159,6 +159,7 @@
#include <cmath>
#include <iosfwd>
#include <iostream> // NOLINT
#include <limits>
#include <string>
#include "Eigen/Core"
@ -197,10 +198,8 @@ struct Jet {
// to be passed in without being fully evaluated until
// they are assigned to v
template<typename Derived>
Jet(const T& value, const Eigen::DenseBase<Derived> &vIn)
: a(value),
v(vIn)
{
EIGEN_STRONG_INLINE Jet(const T& a, const Eigen::DenseBase<Derived> &v)
: a(a), v(v) {
}
// Compound operators
@ -649,7 +648,9 @@ struct NumTraits<ceres::Jet<T, N> > {
return ceres::Jet<T, N>(1e-12);
}
static inline Real epsilon() { return Real(std::numeric_limits<T>::epsilon()); }
static inline Real epsilon() {
return Real(std::numeric_limits<T>::epsilon());
}
enum {
IsComplex = 0,

30
extern/libmv/third_party/ceres/include/ceres/local_parameterization.h vendored
View File

@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
@ -34,6 +34,7 @@
#include <vector>
#include "ceres/internal/port.h"
#include "ceres/internal/disable_warnings.h"
namespace ceres {
@ -109,7 +110,7 @@ namespace ceres {
// Jacobian which is needed to compute the Jacobian of f w.r.t delta.
class CERES_EXPORT LocalParameterization {
public:
virtual ~LocalParameterization() {}
virtual ~LocalParameterization();
// Generalization of the addition operation,
//
@ -121,8 +122,23 @@ class CERES_EXPORT LocalParameterization {
double* x_plus_delta) const = 0;
// The jacobian of Plus(x, delta) w.r.t delta at delta = 0.
//
// jacobian is a row-major GlobalSize() x LocalSize() matrix.
virtual bool ComputeJacobian(const double* x, double* jacobian) const = 0;
// local_matrix = global_matrix * jacobian
//
// global_matrix is a num_rows x GlobalSize row major matrix.
// local_matrix is a num_rows x LocalSize row major matrix.
// jacobian(x) is the matrix returned by ComputeJacobian at x.
//
// This is only used by GradientProblem. For most normal uses, it is
// okay to use the default implementation.
virtual bool MultiplyByJacobian(const double* x,
const int num_rows,
const double* global_matrix,
double* local_matrix) const;
// Size of x.
virtual int GlobalSize() const = 0;
@ -142,6 +158,10 @@ class CERES_EXPORT IdentityParameterization : public LocalParameterization {
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_; }
@ -160,6 +180,10 @@ class CERES_EXPORT SubsetParameterization : public LocalParameterization {
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 static_cast<int>(constancy_mask_.size());
}
@ -188,4 +212,6 @@ class CERES_EXPORT QuaternionParameterization : public LocalParameterization {
} // namespace ceres
#include "ceres/internal/reenable_warnings.h"
#endif // CERES_PUBLIC_LOCAL_PARAMETERIZATION_H_

5
extern/libmv/third_party/ceres/include/ceres/loss_function.h vendored
View File

@ -75,10 +75,11 @@
#ifndef CERES_PUBLIC_LOSS_FUNCTION_H_
#define CERES_PUBLIC_LOSS_FUNCTION_H_
#include "glog/logging.h"
#include "ceres/internal/macros.h"
#include "ceres/internal/scoped_ptr.h"
#include "ceres/types.h"
#include "glog/logging.h"
#include "ceres/internal/disable_warnings.h"
namespace ceres {
@ -395,4 +396,6 @@ class CERES_EXPORT LossFunctionWrapper : public LossFunction {
} // namespace ceres
#include "ceres/internal/disable_warnings.h"
#endif // CERES_PUBLIC_LOSS_FUNCTION_H_

3
extern/libmv/third_party/ceres/include/ceres/normal_prior.h vendored
View File

@ -36,6 +36,7 @@
#include "ceres/cost_function.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/disable_warnings.h"
namespace ceres {
@ -72,4 +73,6 @@ class CERES_EXPORT NormalPrior: public CostFunction {
} // namespace ceres
#include "ceres/internal/reenable_warnings.h"
#endif // CERES_PUBLIC_NORMAL_PRIOR_H_

351
extern/libmv/third_party/ceres/include/ceres/numeric_diff_functor.h vendored
View File

@ -1,351 +0,0 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2013 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
//
// A wrapper class that takes a variadic functor evaluating a
// function, numerically differentiates it and makes it available as a
// templated functor so that it can be easily used as part of Ceres'
// automatic differentiation framework.
//
// For example:
//
// For example, let us assume that
//
// struct IntrinsicProjection
// IntrinsicProjection(const double* observations);
// bool operator()(const double* calibration,
// const double* point,
// double* residuals);
// };
//
// is a functor that implements the projection of a point in its local
// coordinate system onto its image plane and subtracts it from the
// observed point projection.
//
// Now we would like to compose the action of this functor with the
// action of camera extrinsics, i.e., rotation and translation, which
// is given by the following templated function
//
// template<typename T>
// void RotateAndTranslatePoint(const T* rotation,
// const T* translation,
// const T* point,
// T* result);
//
// To compose the extrinsics and intrinsics, we can construct a
// CameraProjection functor as follows.
//
// struct CameraProjection {
// typedef NumericDiffFunctor<IntrinsicProjection, CENTRAL, 2, 5, 3>
// IntrinsicProjectionFunctor;
//
// CameraProjection(double* observation) {
// intrinsic_projection_.reset(
// new IntrinsicProjectionFunctor(observation)) {
// }
//
// template <typename T>
// bool operator()(const T* rotation,
// const T* translation,
// const T* intrinsics,
// const T* point,
// T* residuals) const {
// T transformed_point[3];
// RotateAndTranslatePoint(rotation, translation, point, transformed_point);
// return (*intrinsic_projection_)(intrinsics, transformed_point, residual);
// }
//
// private:
// scoped_ptr<IntrinsicProjectionFunctor> intrinsic_projection_;
// };
//
// Here, we made the choice of using CENTRAL differences to compute
// the jacobian of IntrinsicProjection.
//
// Now, we are ready to construct an automatically differentiated cost
// function as
//
// CostFunction* cost_function =
// new AutoDiffCostFunction<CameraProjection, 2, 3, 3, 5>(
// new CameraProjection(observations));
//
// cost_function now seamlessly integrates automatic differentiation
// of RotateAndTranslatePoint with a numerically differentiated
// version of IntrinsicProjection.
#ifndef CERES_PUBLIC_NUMERIC_DIFF_FUNCTOR_H_
#define CERES_PUBLIC_NUMERIC_DIFF_FUNCTOR_H_
#include "ceres/numeric_diff_cost_function.h"
#include "ceres/types.h"
#include "ceres/cost_function_to_functor.h"
namespace ceres {
template<typename Functor,
NumericDiffMethod kMethod = CENTRAL,
int kNumResiduals = 0,
int N0 = 0, int N1 = 0 , int N2 = 0, int N3 = 0, int N4 = 0,
int N5 = 0, int N6 = 0 , int N7 = 0, int N8 = 0, int N9 = 0>
class NumericDiffFunctor {
public:
// relative_step_size controls the step size used by the numeric
// differentiation process.
explicit NumericDiffFunctor(double relative_step_size = 1e-6)
: functor_(
new NumericDiffCostFunction<Functor,
kMethod,
kNumResiduals,
N0, N1, N2, N3, N4,
N5, N6, N7, N8, N9>(new Functor,
TAKE_OWNERSHIP,
kNumResiduals,
relative_step_size)) {
}
NumericDiffFunctor(Functor* functor, double relative_step_size = 1e-6)
: functor_(new NumericDiffCostFunction<Functor,
kMethod,
kNumResiduals,
N0, N1, N2, N3, N4,
N5, N6, N7, N8, N9>(
functor,
TAKE_OWNERSHIP,
kNumResiduals,
relative_step_size)) {
}
bool operator()(const double* x0, double* residuals) const {
return functor_(x0, residuals);
}
bool operator()(const double* x0,
const double* x1,
double* residuals) const {
return functor_(x0, x1, residuals);
}
bool operator()(const double* x0,
const double* x1,
const double* x2,
double* residuals) const {
return functor_(x0, x1, x2, residuals);
}
bool operator()(const double* x0,
const double* x1,
const double* x2,
const double* x3,
double* residuals) const {
return functor_(x0, x1, x2, x3, residuals);
}
bool operator()(const double* x0,
const double* x1,
const double* x2,
const double* x3,
const double* x4,
double* residuals) const {
return functor_(x0, x1, x2, x3, x4, residuals);
}
bool operator()(const double* x0,
const double* x1,
const double* x2,
const double* x3,
const double* x4,
const double* x5,
double* residuals) const {
return functor_(x0, x1, x2, x3, x4, x5, residuals);
}
bool operator()(const double* x0,
const double* x1,
const double* x2,
const double* x3,
const double* x4,
const double* x5,
const double* x6,
double* residuals) const {
return functor_(x0, x1, x2, x3, x4, x5, x6, residuals);
}
bool operator()(const double* x0,
const double* x1,
const double* x2,
const double* x3,
const double* x4,
const double* x5,
const double* x6,
const double* x7,
double* residuals) const {
return functor_(x0, x1, x2, x3, x4, x5, x6, x7, residuals);
}
bool operator()(const double* x0,
const double* x1,
const double* x2,
const double* x3,
const double* x4,
const double* x5,
const double* x6,
const double* x7,
const double* x8,
double* residuals) const {
return functor_(x0, x1, x2, x3, x4, x5, x6, x7, x8, residuals);
}
bool operator()(const double* x0,
const double* x1,
const double* x2,
const double* x3,
const double* x4,
const double* x5,
const double* x6,
const double* x7,
const double* x8,
const double* x9,
double* residuals) const {
return functor_(x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, residuals);
}
template <typename T>
bool operator()(const T* x0, T* residuals) const {
return functor_(x0, residuals);
}
template <typename T>
bool operator()(const T* x0,
const T* x1,
T* residuals) const {
return functor_(x0, x1, residuals);
}
template <typename T>
bool operator()(const T* x0,
const T* x1,
const T* x2,
T* residuals) const {
return functor_(x0, x1, x2, residuals);
}
template <typename T>
bool operator()(const T* x0,
const T* x1,
const T* x2,
const T* x3,
T* residuals) const {
return functor_(x0, x1, x2, x3, residuals);
}
template <typename T>
bool operator()(const T* x0,
const T* x1,
const T* x2,
const T* x3,
const T* x4,
T* residuals) const {
return functor_(x0, x1, x2, x3, x4, residuals);
}
template <typename T>
bool operator()(const T* x0,
const T* x1,
const T* x2,
const T* x3,
const T* x4,
const T* x5,
T* residuals) const {
return functor_(x0, x1, x2, x3, x4, x5, residuals);
}
template <typename T>
bool operator()(const T* x0,
const T* x1,
const T* x2,
const T* x3,
const T* x4,
const T* x5,
const T* x6,
T* residuals) const {
return functor_(x0, x1, x2, x3, x4, x5, x6, residuals);
}
template <typename T>
bool operator()(const T* x0,
const T* x1,
const T* x2,
const T* x3,
const T* x4,
const T* x5,
const T* x6,
const T* x7,
T* residuals) const {
return functor_(x0, x1, x2, x3, x4, x5, x6, x7, residuals);
}
template <typename T>
bool operator()(const T* x0,
const T* x1,
const T* x2,
const T* x3,
const T* x4,
const T* x5,
const T* x6,
const T* x7,
const T* x8,
T* residuals) const {
return functor_(x0, x1, x2, x3, x4, x5, x6, x7, x8, residuals);
}
template <typename T>
bool operator()(const T* x0,
const T* x1,
const T* x2,
const T* x3,
const T* x4,
const T* x5,
const T* x6,
const T* x7,
const T* x8,
const T* x9,
T* residuals) const {
return functor_(x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, residuals);
}
private:
CostFunctionToFunctor<kNumResiduals,
N0, N1, N2, N3, N4,
N5, N6, N7, N8, N9> functor_;
};
} // namespace ceres
#endif // CERES_PUBLIC_NUMERIC_DIFF_FUNCTOR_H_

28
extern/libmv/third_party/ceres/include/ceres/ordered_groups.h vendored
View File

@ -33,7 +33,9 @@
#include <map>
#include <set>
#include <vector>
#include "ceres/internal/port.h"
#include "glog/logging.h"
namespace ceres {
@ -103,6 +105,20 @@ class OrderedGroups {
return true;
}
// Bulk remove elements. The return value indicates the number of
// elements successfully removed.
int Remove(const vector<T>& elements) {
if (NumElements() == 0 || elements.size() == 0) {
return 0;
}
int num_removed = 0;
for (int i = 0; i < elements.size(); ++i) {
num_removed += Remove(elements[i]);
}
return num_removed;
}
// Reverse the order of the groups in place.
void Reverse() {
typename map<int, set<T> >::reverse_iterator it =
@ -156,10 +172,22 @@ class OrderedGroups {
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
// crash.
int MinNonZeroGroup() const {
CHECK_NE(NumGroups(), 0);
return group_to_elements_.begin()->first;
}
const map<int, set<T> >& group_to_elements() const {
return group_to_elements_;
}
const map<T, int>& element_to_group() const {
return element_to_group_;
}
private:
map<int, set<T> > group_to_elements_;
map<T, int> element_to_group_;

14
extern/libmv/third_party/ceres/include/ceres/problem.h vendored
View File

@ -39,11 +39,12 @@
#include <set>
#include <vector>
#include "glog/logging.h"
#include "ceres/internal/macros.h"
#include "ceres/internal/port.h"
#include "ceres/internal/scoped_ptr.h"
#include "ceres/types.h"
#include "glog/logging.h"
#include "ceres/internal/disable_warnings.h"
namespace ceres {
@ -367,6 +368,15 @@ class CERES_EXPORT Problem {
const ResidualBlockId residual_block,
vector<double*>* parameter_blocks) const;
// Get the CostFunction for the given residual block.
const CostFunction* GetCostFunctionForResidualBlock(
const ResidualBlockId residual_block) const;
// Get the LossFunction for the given residual block. Returns NULL
// if no loss function is associated with this residual block.
const LossFunction* GetLossFunctionForResidualBlock(
const ResidualBlockId residual_block) const;
// Get all the residual blocks that depend on the given parameter block.
//
// If Problem::Options::enable_fast_removal is true, then
@ -466,4 +476,6 @@ class CERES_EXPORT Problem {
} // namespace ceres
#include "ceres/internal/reenable_warnings.h"
#endif // CERES_PUBLIC_PROBLEM_H_

48
extern/libmv/third_party/ceres/include/ceres/solver.h vendored
View File

@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
@ -40,6 +40,7 @@
#include "ceres/iteration_callback.h"
#include "ceres/ordered_groups.h"
#include "ceres/types.h"
#include "ceres/internal/disable_warnings.h"
namespace ceres {
@ -91,7 +92,7 @@ class CERES_EXPORT Solver {
gradient_tolerance = 1e-10;
parameter_tolerance = 1e-8;
#if defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CXSPARSE)
#if defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CXSPARSE) && !defined(CERES_ENABLE_LGPL_CODE)
linear_solver_type = DENSE_QR;
#else
linear_solver_type = SPARSE_NORMAL_CHOLESKY;
@ -100,16 +101,26 @@ class CERES_EXPORT Solver {
preconditioner_type = JACOBI;
visibility_clustering_type = CANONICAL_VIEWS;
dense_linear_algebra_library_type = EIGEN;
sparse_linear_algebra_library_type = SUITE_SPARSE;
#if defined(CERES_NO_SUITESPARSE) && !defined(CERES_NO_CXSPARSE)
sparse_linear_algebra_library_type = CX_SPARSE;
#endif
// Choose a default sparse linear algebra library in the order:
//
// SUITE_SPARSE > CX_SPARSE > EIGEN_SPARSE
#if !defined(CERES_NO_SUITESPARSE)
sparse_linear_algebra_library_type = SUITE_SPARSE;
#else
#if !defined(CERES_NO_CXSPARSE)
sparse_linear_algebra_library_type = CX_SPARSE;
#else
#if defined(CERES_USE_EIGEN_SPARSE)
sparse_linear_algebra_library_type = EIGEN_SPARSE;
#endif
#endif
#endif
num_linear_solver_threads = 1;
use_postordering = false;
dynamic_sparsity = false;
min_linear_solver_iterations = 1;
min_linear_solver_iterations = 0;
max_linear_solver_iterations = 500;
eta = 1e-1;
jacobi_scaling = true;
@ -125,6 +136,11 @@ class CERES_EXPORT Solver {
update_state_every_iteration = false;
}
// Returns true if the options struct has a valid
// configuration. Returns false otherwise, and fills in *error
// with a message describing the problem.
bool IsValid(string* error) const;
// Minimizer options ----------------------------------------
// Ceres supports the two major families of optimization strategies -
@ -707,10 +723,6 @@ class CERES_EXPORT Solver {
//
// The solver does NOT take ownership of these pointers.
vector<IterationCallback*> callbacks;
// If non-empty, a summary of the execution of the solver is
// recorded to this file.
string solver_log;
};
struct CERES_EXPORT Summary {
@ -898,9 +910,15 @@ class CERES_EXPORT Solver {
// parameter blocks.
vector<int> inner_iteration_ordering_used;
// Type of preconditioner used for solving the trust region
// step. Only meaningful when an iterative linear solver is used.
PreconditionerType preconditioner_type;
// Type of the preconditioner requested by the user.
PreconditionerType preconditioner_type_given;
// Type of the preconditioner actually used. This may be different
// from linear_solver_type_given if Ceres determines that the
// problem structure is not compatible with the linear solver
// requested or if the linear solver requested by the user is not
// available.
PreconditionerType preconditioner_type_used;
// Type of clustering algorithm used for visibility based
// preconditioning. Only meaningful when the preconditioner_type
@ -957,4 +975,6 @@ CERES_EXPORT void Solve(const Solver::Options& options,
} // namespace ceres
#include "ceres/internal/reenable_warnings.h"
#endif // CERES_PUBLIC_SOLVER_H_

19
extern/libmv/third_party/ceres/include/ceres/types.h vendored
View File

@ -40,6 +40,7 @@
#include <string>
#include "ceres/internal/port.h"
#include "ceres/internal/disable_warnings.h"
namespace ceres {
@ -149,8 +150,14 @@ enum SparseLinearAlgebraLibraryType {
// minimum degree ordering.
SUITE_SPARSE,
// A lightweight replacment for SuiteSparse.
CX_SPARSE
// A lightweight replacment for SuiteSparse, which does not require
// a LAPACK/BLAS implementation. Consequently, its performance is
// also a bit lower than SuiteSparse.
CX_SPARSE,
// Eigen's sparse linear algebra routines. In particular Ceres uses
// the Simplicial LDLT routines.
EIGEN_SPARSE
};
enum DenseLinearAlgebraLibraryType {
@ -246,7 +253,7 @@ enum LineSearchDirectionType {
// details see Numerical Optimization by Nocedal & Wright.
enum NonlinearConjugateGradientType {
FLETCHER_REEVES,
POLAK_RIBIRERE,
POLAK_RIBIERE,
HESTENES_STIEFEL,
};
@ -398,8 +405,8 @@ enum LineSearchInterpolationType {
enum CovarianceAlgorithmType {
DENSE_SVD,
SPARSE_CHOLESKY,
SPARSE_QR
SUITE_SPARSE_QR,
EIGEN_SPARSE_QR
};
CERES_EXPORT const char* LinearSolverTypeToString(
@ -475,4 +482,6 @@ CERES_EXPORT bool IsDenseLinearAlgebraLibraryTypeAvailable(
} // namespace ceres
#include "ceres/internal/reenable_warnings.h"
#endif // CERES_PUBLIC_TYPES_H_

49
extern/libmv/third_party/ceres/include/ceres/version.h vendored Normal file
View File

@ -0,0 +1,49 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 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: mierle@gmail.com (Keir Mierle)
#ifndef CERES_PUBLIC_VERSION_H_
#define CERES_PUBLIC_VERSION_H_
#define CERES_VERSION_MAJOR 1
#define CERES_VERSION_MINOR 10
#define CERES_VERSION_REVISION 0
#define CERES_VERSION_ABI 1
// Classic CPP stringifcation; the extra level of indirection allows the
// preprocessor to expand the macro before being converted to a string.
#define CERES_TO_STRING_HELPER(x) #x
#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)
#endif // CERES_PUBLIC_VERSION_H_

287
extern/libmv/third_party/ceres/internal/ceres/CMakeLists.txt vendored
View File

@ -1,287 +0,0 @@
# Ceres Solver - A fast non-linear least squares minimizer
# Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
# http://code.google.com/p/ceres-solver/
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
# * Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimer in the documentation
# and/or other materials provided with the distribution.
# * Neither the name of Google Inc. nor the names of its contributors may be
# used to endorse or promote products derived from this software without
# specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
#
# Author: keir@google.com (Keir Mierle)
SET(CERES_INTERNAL_SRC
array_utils.cc
blas.cc
block_evaluate_preparer.cc
block_jacobi_preconditioner.cc
block_jacobian_writer.cc
block_random_access_dense_matrix.cc
block_random_access_diagonal_matrix.cc
block_random_access_matrix.cc
block_random_access_sparse_matrix.cc
block_sparse_matrix.cc
block_structure.cc
c_api.cc
canonical_views_clustering.cc
cgnr_solver.cc
compressed_col_sparse_matrix_utils.cc
compressed_row_jacobian_writer.cc
compressed_row_sparse_matrix.cc
conditioned_cost_function.cc
conjugate_gradients_solver.cc
coordinate_descent_minimizer.cc
corrector.cc
covariance.cc
covariance_impl.cc
cxsparse.cc
dense_normal_cholesky_solver.cc
dense_qr_solver.cc
dense_sparse_matrix.cc
detect_structure.cc
dogleg_strategy.cc
dynamic_compressed_row_jacobian_writer.cc
dynamic_compressed_row_sparse_matrix.cc
evaluator.cc
file.cc
gradient_checking_cost_function.cc
implicit_schur_complement.cc
incomplete_lq_factorization.cc
iterative_schur_complement_solver.cc
levenberg_marquardt_strategy.cc
lapack.cc
line_search.cc
line_search_direction.cc
line_search_minimizer.cc
linear_least_squares_problems.cc
linear_operator.cc
linear_solver.cc
local_parameterization.cc
loss_function.cc
low_rank_inverse_hessian.cc
minimizer.cc
normal_prior.cc
parameter_block_ordering.cc
partitioned_matrix_view.cc
polynomial.cc
preconditioner.cc
problem.cc
problem_impl.cc
program.cc
residual_block.cc
residual_block_utils.cc
schur_complement_solver.cc
schur_eliminator.cc
schur_jacobi_preconditioner.cc
scratch_evaluate_preparer.cc
single_linkage_clustering.cc
solver.cc
solver_impl.cc
sparse_matrix.cc
sparse_normal_cholesky_solver.cc
split.cc
stringprintf.cc
suitesparse.cc
triplet_sparse_matrix.cc
trust_region_minimizer.cc
trust_region_strategy.cc
types.cc
visibility.cc
visibility_based_preconditioner.cc
wall_time.cc
)
# Heuristic for determining LIB_SUFFIX. FHS recommends that 64-bit systems
# install native libraries to lib64 rather than lib. Most distros seem to
# follow this convention with a couple notable exceptions (Debian-based and
# Arch-based distros) which we try to detect here.
IF (CMAKE_SYSTEM_NAME MATCHES "Linux" AND
NOT DEFINED LIB_SUFFIX AND
NOT CMAKE_CROSSCOMPILING AND
CMAKE_SIZEOF_VOID_P EQUAL "8" AND
NOT EXISTS "/etc/debian_version" AND
NOT EXISTS "/etc/arch-release")
SET(LIB_SUFFIX "64")
ENDIF ()
# Also depend on the header files so that they appear in IDEs.
FILE(GLOB CERES_INTERNAL_HDRS *.h)
# Include the specialized schur solvers.
IF (SCHUR_SPECIALIZATIONS)
FILE(GLOB CERES_INTERNAL_SCHUR_FILES generated/*.cc)
ELSE (SCHUR_SPECIALIZATIONS)
# Only the fully dynamic solver. The build is much faster this way.
FILE(GLOB CERES_INTERNAL_SCHUR_FILES generated/*_d_d_d.cc)
ENDIF (SCHUR_SPECIALIZATIONS)
# Primarily for Android, but optionally for others, use the minimal internal
# Glog implementation.
IF (MINIGLOG)
ADD_LIBRARY(miniglog STATIC miniglog/glog/logging.cc)
INSTALL(TARGETS miniglog
EXPORT CeresExport
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib${LIB_SUFFIX}
ARCHIVE DESTINATION lib${LIB_SUFFIX})
ENDIF (MINIGLOG)
SET(CERES_LIBRARY_PUBLIC_DEPENDENCIES ${GLOG_LIBRARIES})
IF (SUITESPARSE AND SUITESPARSE_FOUND)
LIST(APPEND CERES_LIBRARY_PRIVATE_DEPENDENCIES ${SUITESPARSE_LIBRARIES})
ENDIF (SUITESPARSE AND SUITESPARSE_FOUND)
IF (CXSPARSE AND CXSPARSE_FOUND)
LIST(APPEND CERES_LIBRARY_PRIVATE_DEPENDENCIES ${CXSPARSE_LIBRARIES})
ENDIF (CXSPARSE AND CXSPARSE_FOUND)
IF (BLAS_FOUND AND LAPACK_FOUND)
LIST(APPEND CERES_LIBRARY_PRIVATE_DEPENDENCIES ${LAPACK_LIBRARIES})
LIST(APPEND CERES_LIBRARY_PRIVATE_DEPENDENCIES ${BLAS_LIBRARIES})
ENDIF (BLAS_FOUND AND LAPACK_FOUND)
IF (OPENMP_FOUND)
IF (NOT MSVC)
LIST(APPEND CERES_LIBRARY_PRIVATE_DEPENDENCIES gomp)
LIST(APPEND CERES_LIBRARY_PRIVATE_DEPENDENCIES ${CMAKE_THREAD_LIBS_INIT})
ENDIF (NOT MSVC)
ENDIF (OPENMP_FOUND)
SET(CERES_LIBRARY_SOURCE
${CERES_INTERNAL_SRC}
${CERES_INTERNAL_HDRS}
${CERES_INTERNAL_SCHUR_FILES})
ADD_LIBRARY(ceres ${CERES_LIBRARY_SOURCE})
SET_TARGET_PROPERTIES(ceres PROPERTIES
VERSION ${CERES_VERSION}
SOVERSION ${CERES_VERSION_MAJOR}
)
IF (BUILD_SHARED_LIBS)
# When building a shared library, mark all external libraries as
# PRIVATE so they don't show up as a dependency.
TARGET_LINK_LIBRARIES(ceres
LINK_PUBLIC ${CERES_LIBRARY_PUBLIC_DEPENDENCIES}
LINK_PRIVATE ${CERES_LIBRARY_PRIVATE_DEPENDENCIES})
ELSE (BUILD_SHARED_LIBS)
# When building a static library, all external libraries are
# PUBLIC(default) since the user needs to link to them.
# They will be listed in CeresTargets.cmake.
SET(CERES_LIBRARY_DEPENDENCIES
${CERES_LIBRARY_PUBLIC_DEPENDENCIES}
${CERES_LIBRARY_PRIVATE_DEPENDENCIES})
TARGET_LINK_LIBRARIES(ceres ${CERES_LIBRARY_DEPENDENCIES})
ENDIF (BUILD_SHARED_LIBS)
INSTALL(TARGETS ceres
EXPORT CeresExport
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib${LIB_SUFFIX}
ARCHIVE DESTINATION lib${LIB_SUFFIX})
IF (BUILD_TESTING AND GFLAGS)
ADD_LIBRARY(gtest gmock_gtest_all.cc gmock_main.cc)
ADD_LIBRARY(test_util
evaluator_test_utils.cc
numeric_diff_test_utils.cc
test_util.cc)
TARGET_LINK_LIBRARIES(gtest ${GFLAGS_LIBRARIES} ${GLOG_LIBRARIES})
TARGET_LINK_LIBRARIES(test_util ceres gtest ${GLOG_LIBRARIES})
MACRO (CERES_TEST NAME)
ADD_EXECUTABLE(${NAME}_test ${NAME}_test.cc)
TARGET_LINK_LIBRARIES(${NAME}_test test_util ceres gtest)
ADD_TEST(NAME ${NAME}_test
COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/${NAME}_test
--test_srcdir
${CMAKE_SOURCE_DIR}/data)
ENDMACRO (CERES_TEST)
CERES_TEST(array_utils)
CERES_TEST(autodiff)
CERES_TEST(autodiff_cost_function)
CERES_TEST(autodiff_local_parameterization)
CERES_TEST(block_random_access_dense_matrix)
CERES_TEST(block_random_access_diagonal_matrix)
CERES_TEST(block_random_access_sparse_matrix)
CERES_TEST(block_sparse_matrix)
CERES_TEST(c_api)
CERES_TEST(canonical_views_clustering)
CERES_TEST(compressed_row_sparse_matrix)
CERES_TEST(conditioned_cost_function)
CERES_TEST(corrector)
CERES_TEST(cost_function_to_functor)
CERES_TEST(covariance)
CERES_TEST(dense_sparse_matrix)
CERES_TEST(dynamic_autodiff_cost_function)
CERES_TEST(dynamic_compressed_row_sparse_matrix)
CERES_TEST(dynamic_numeric_diff_cost_function)
CERES_TEST(evaluator)
CERES_TEST(gradient_checker)
CERES_TEST(gradient_checking_cost_function)
CERES_TEST(graph)
CERES_TEST(graph_algorithms)
CERES_TEST(implicit_schur_complement)
CERES_TEST(incomplete_lq_factorization)
CERES_TEST(iterative_schur_complement_solver)
CERES_TEST(jet)
CERES_TEST(levenberg_marquardt_strategy)
CERES_TEST(dogleg_strategy)
CERES_TEST(local_parameterization)
CERES_TEST(loss_function)
CERES_TEST(minimizer)
CERES_TEST(normal_prior)
CERES_TEST(numeric_diff_cost_function)
CERES_TEST(numeric_diff_functor)
CERES_TEST(ordered_groups)
CERES_TEST(parameter_block)
CERES_TEST(parameter_block_ordering)
CERES_TEST(partitioned_matrix_view)
CERES_TEST(polynomial)
CERES_TEST(problem)
CERES_TEST(residual_block)
CERES_TEST(residual_block_utils)
CERES_TEST(rotation)
CERES_TEST(schur_complement_solver)
CERES_TEST(schur_eliminator)
CERES_TEST(single_linkage_clustering)
CERES_TEST(small_blas)
CERES_TEST(solver_impl)
# TODO(sameeragarwal): This test should ultimately be made
# independent of SuiteSparse.
IF (SUITESPARSE AND SUITESPARSE_FOUND)
CERES_TEST(compressed_col_sparse_matrix_utils)
ENDIF (SUITESPARSE AND SUITESPARSE_FOUND)
CERES_TEST(symmetric_linear_solver)
CERES_TEST(triplet_sparse_matrix)
CERES_TEST(trust_region_minimizer)
CERES_TEST(unsymmetric_linear_solver)
CERES_TEST(visibility)
CERES_TEST(visibility_based_preconditioner)
# Put the large end to end test last.
CERES_TEST(system)
ENDIF (BUILD_TESTING AND GFLAGS)

17
extern/libmv/third_party/ceres/internal/ceres/array_utils.cc vendored
View File

@ -30,10 +30,11 @@
#include "ceres/array_utils.h"
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <string>
#include <vector>
#include "ceres/fpclassify.h"
#include "ceres/stringprintf.h"
@ -94,5 +95,19 @@ void AppendArrayToString(const int size, const double* x, string* result) {
}
}
void MapValuesToContiguousRange(const int size, int* array) {
std::vector<int> unique_values(array, array + size);
std::sort(unique_values.begin(), unique_values.end());
unique_values.erase(std::unique(unique_values.begin(),
unique_values.end()),
unique_values.end());
for (int i = 0; i < size; ++i) {
array[i] = std::lower_bound(unique_values.begin(),
unique_values.end(),
array[i]) - unique_values.begin();
}
}
} // namespace internal
} // namespace ceres

15
extern/libmv/third_party/ceres/internal/ceres/array_utils.h vendored
View File

@ -67,6 +67,21 @@ void AppendArrayToString(const int size, const double* x, string* result);
extern const double kImpossibleValue;
// This routine takes an array of integer values, sorts and uniques
// them and then maps each value in the array to its position in the
// sorted+uniqued array. By doing this, if there are are k unique
// values in the array, each value is replaced by an integer in the
// range [0, k-1], while preserving their relative order.
//
// For example
//
// [1 0 3 5 0 1 5]
//
// gets mapped to
//
// [1 0 2 3 0 1 3]
void MapValuesToContiguousRange(int size, int* array);
} // namespace internal
} // namespace ceres

109
extern/libmv/third_party/ceres/internal/ceres/callbacks.cc vendored Normal file
View File

@ -0,0 +1,109 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#include <iostream> // NO LINT
#include "ceres/callbacks.h"
#include "ceres/program.h"
#include "ceres/stringprintf.h"
#include "glog/logging.h"
namespace ceres {
namespace internal {
StateUpdatingCallback::StateUpdatingCallback(Program* program,
double* parameters)
: program_(program), parameters_(parameters) {}
StateUpdatingCallback::~StateUpdatingCallback() {}
CallbackReturnType StateUpdatingCallback::operator()(
const IterationSummary& summary) {
if (summary.step_is_successful) {
program_->StateVectorToParameterBlocks(parameters_);
program_->CopyParameterBlockStateToUserState();
}
return SOLVER_CONTINUE;
}
LoggingCallback::LoggingCallback(const MinimizerType minimizer_type,
const bool log_to_stdout)
: minimizer_type(minimizer_type),
log_to_stdout_(log_to_stdout) {}
LoggingCallback::~LoggingCallback() {}
CallbackReturnType LoggingCallback::operator()(
const IterationSummary& summary) {
string output;
if (minimizer_type == LINE_SEARCH) {
const char* kReportRowFormat =
"% 4d: f:% 8e d:% 3.2e g:% 3.2e h:% 3.2e "
"s:% 3.2e e:% 3d it:% 3.2e tt:% 3.2e";
output = StringPrintf(kReportRowFormat,
summary.iteration,
summary.cost,
summary.cost_change,
summary.gradient_max_norm,
summary.step_norm,
summary.step_size,
summary.line_search_function_evaluations,
summary.iteration_time_in_seconds,
summary.cumulative_time_in_seconds);
} else if (minimizer_type == TRUST_REGION) {
if (summary.iteration == 0) {
output = "iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time\n";
}
const char* kReportRowFormat =
"% 4d % 8e % 3.2e % 3.2e % 3.2e % 3.2e % 3.2e % 3d % 3.2e % 3.2e";
output += StringPrintf(kReportRowFormat,
summary.iteration,
summary.cost,
summary.cost_change,
summary.gradient_max_norm,
summary.step_norm,
summary.relative_decrease,
summary.trust_region_radius,
summary.linear_solver_iterations,
summary.iteration_time_in_seconds,
summary.cumulative_time_in_seconds);
} else {
LOG(FATAL) << "Unknown minimizer type.";
}
if (log_to_stdout_) {
cout << output << endl;
} else {
VLOG(1) << output;
}
return SOLVER_CONTINUE;
}
} // namespace internal
} // namespace ceres

71
extern/libmv/third_party/ceres/internal/ceres/callbacks.h vendored Normal file
View File

@ -0,0 +1,71 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#ifndef CERES_INTERNAL_CALLBACKS_H_
#define CERES_INTERNAL_CALLBACKS_H_
#include <string>
#include "ceres/iteration_callback.h"
#include "ceres/internal/port.h"
namespace ceres {
namespace internal {
class Program;
// Callback for updating the externally visible state of parameter
// blocks.
class StateUpdatingCallback : public IterationCallback {
public:
StateUpdatingCallback(Program* program, double* parameters);
virtual ~StateUpdatingCallback();
virtual CallbackReturnType operator()(const IterationSummary& summary);
private:
Program* program_;
double* parameters_;
};
// Callback for logging the state of the minimizer to STDERR or
// STDOUT depending on the user's preferences and logging level.
class LoggingCallback : public IterationCallback {
public:
LoggingCallback(MinimizerType minimizer_type, bool log_to_stdout);
virtual ~LoggingCallback();
virtual CallbackReturnType operator()(const IterationSummary& summary);
private:
const MinimizerType minimizer_type;
const bool log_to_stdout_;
};
} // namespace internal
} // namespace ceres
#endif // CERES_INTERNAL_CALLBACKS_H_

10
extern/libmv/third_party/ceres/internal/ceres/canonical_views_clustering.cc vendored
View File

@ -61,7 +61,7 @@ class CanonicalViewsClustering {
// vertices may not be assigned to any cluster. In this case they
// are assigned to a cluster with id = kInvalidClusterId.
void ComputeClustering(const CanonicalViewsClusteringOptions& options,
const Graph<int>& graph,
const WeightedGraph<int>& graph,
vector<int>* centers,
IntMap* membership);
@ -74,7 +74,7 @@ class CanonicalViewsClustering {
IntMap* membership) const;
CanonicalViewsClusteringOptions options_;
const Graph<int>* graph_;
const WeightedGraph<int>* graph_;
// Maps a view to its representative canonical view (its cluster
// center).
IntMap view_to_canonical_view_;
@ -85,7 +85,7 @@ class CanonicalViewsClustering {
void ComputeCanonicalViewsClustering(
const CanonicalViewsClusteringOptions& options,
const Graph<int>& graph,
const WeightedGraph<int>& graph,
vector<int>* centers,
IntMap* membership) {
time_t start_time = time(NULL);
@ -98,7 +98,7 @@ void ComputeCanonicalViewsClustering(
// Implementation of CanonicalViewsClustering
void CanonicalViewsClustering::ComputeClustering(
const CanonicalViewsClusteringOptions& options,
const Graph<int>& graph,
const WeightedGraph<int>& graph,
vector<int>* centers,
IntMap* membership) {
options_ = options;
@ -150,7 +150,7 @@ void CanonicalViewsClustering::FindValidViews(
for (IntSet::const_iterator view = views.begin();
view != views.end();
++view) {
if (graph_->VertexWeight(*view) != Graph<int>::InvalidWeight()) {
if (graph_->VertexWeight(*view) != WeightedGraph<int>::InvalidWeight()) {
valid_views->insert(*view);
}
}

2
extern/libmv/third_party/ceres/internal/ceres/canonical_views_clustering.h vendored
View File

@ -101,7 +101,7 @@ struct CanonicalViewsClusteringOptions;
// cluster. In this case they are assigned to a cluster with id = -1;
void ComputeCanonicalViewsClustering(
const CanonicalViewsClusteringOptions& options,
const Graph<int>& graph,
const WeightedGraph<int>& graph,
vector<int>* centers,
HashMap<int, int>* membership);

17
extern/libmv/third_party/ceres/internal/ceres/conjugate_gradients_solver.cc vendored
View File

@ -101,7 +101,7 @@ LinearSolver::Summary ConjugateGradientsSolver::Solve(
A->RightMultiply(x, tmp.data());
r = bref - tmp;
double norm_r = r.norm();
if (norm_r <= tol_r) {
if (options_.min_num_iterations == 0 && norm_r <= tol_r) {
summary.termination_type = LINEAR_SOLVER_SUCCESS;
summary.message =
StringPrintf("Convergence. |r| = %e <= %e.", norm_r, tol_r);
@ -113,9 +113,8 @@ LinearSolver::Summary ConjugateGradientsSolver::Solve(
// Initial value of the quadratic model Q = x'Ax - 2 * b'x.
double Q0 = -1.0 * xref.dot(bref + r);
for (summary.num_iterations = 1;
summary.num_iterations < options_.max_num_iterations;
++summary.num_iterations) {
for (summary.num_iterations = 1;; ++summary.num_iterations) {
// Apply preconditioner
if (per_solve_options.preconditioner != NULL) {
z.setZero();
@ -207,7 +206,8 @@ LinearSolver::Summary ConjugateGradientsSolver::Solve(
// 124(1-2), 45-59, 2000.
//
const double zeta = summary.num_iterations * (Q1 - Q0) / Q1;
if (zeta < per_solve_options.q_tolerance) {
if (zeta < per_solve_options.q_tolerance &&
summary.num_iterations >= options_.min_num_iterations) {
summary.termination_type = LINEAR_SOLVER_SUCCESS;
summary.message =
StringPrintf("Convergence: zeta = %e < %e",
@ -219,12 +219,17 @@ LinearSolver::Summary ConjugateGradientsSolver::Solve(
// Residual based termination.
norm_r = r. norm();
if (norm_r <= tol_r) {
if (norm_r <= tol_r &&
summary.num_iterations >= options_.min_num_iterations) {
summary.termination_type = LINEAR_SOLVER_SUCCESS;
summary.message =
StringPrintf("Convergence. |r| = %e <= %e.", norm_r, tol_r);
break;
}
if (summary.num_iterations >= options_.max_num_iterations) {
break;
}
}
return summary;

56
extern/libmv/third_party/ceres/internal/ceres/coordinate_descent_minimizer.cc vendored
View File

@ -40,15 +40,15 @@
#include "ceres/evaluator.h"
#include "ceres/linear_solver.h"
#include "ceres/minimizer.h"
#include "ceres/ordered_groups.h"
#include "ceres/parameter_block.h"
#include "ceres/parameter_block_ordering.h"
#include "ceres/problem_impl.h"
#include "ceres/program.h"
#include "ceres/residual_block.h"
#include "ceres/solver.h"
#include "ceres/solver_impl.h"
#include "ceres/trust_region_minimizer.h"
#include "ceres/trust_region_strategy.h"
#include "ceres/parameter_block_ordering.h"
namespace ceres {
namespace internal {
@ -210,28 +210,54 @@ void CoordinateDescentMinimizer::Solve(Program* program,
summary->final_cost = 0.0;
string error;
scoped_ptr<Evaluator> evaluator(
Evaluator::Create(evaluator_options_, program, &error));
CHECK_NOTNULL(evaluator.get());
scoped_ptr<SparseMatrix> jacobian(evaluator->CreateJacobian());
CHECK_NOTNULL(jacobian.get());
Minimizer::Options minimizer_options;
minimizer_options.evaluator.reset(
CHECK_NOTNULL(Evaluator::Create(evaluator_options_, program, &error)));
minimizer_options.jacobian.reset(
CHECK_NOTNULL(minimizer_options.evaluator->CreateJacobian()));
TrustRegionStrategy::Options trs_options;
trs_options.linear_solver = linear_solver;
scoped_ptr<TrustRegionStrategy>trust_region_strategy(
minimizer_options.trust_region_strategy.reset(
CHECK_NOTNULL(TrustRegionStrategy::Create(trs_options)));
Minimizer::Options minimizer_options;
minimizer_options.evaluator = evaluator.get();
minimizer_options.jacobian = jacobian.get();
minimizer_options.trust_region_strategy = trust_region_strategy.get();
minimizer_options.is_silent = true;
TrustRegionMinimizer minimizer;
minimizer.Minimize(minimizer_options, parameter, summary);
}
bool CoordinateDescentMinimizer::IsOrderingValid(
const Program& program,
const ParameterBlockOrdering& ordering,
string* message) {
const map<int, set<double*> >& group_to_elements =
ordering.group_to_elements();
// Verify that each group is an independent set
map<int, set<double*> >::const_iterator it = group_to_elements.begin();
for ( ; it != group_to_elements.end(); ++it) {
if (!program.IsParameterBlockSetIndependent(it->second)) {
*message =
StringPrintf("The user-provided "
"parameter_blocks_for_inner_iterations does not "
"form an independent set. Group Id: %d", it->first);
return false;
}
}
return true;
}
// Find a recursive decomposition of the Hessian matrix as a set
// of independent sets of decreasing size and invert it. This
// seems to work better in practice, i.e., Cameras before
// points.
ParameterBlockOrdering* CoordinateDescentMinimizer::CreateOrdering(
const Program& program) {
scoped_ptr<ParameterBlockOrdering> ordering(new ParameterBlockOrdering);
ComputeRecursiveIndependentSetOrdering(program, ordering.get());
ordering->Reverse();
return ordering.release();
}
} // namespace internal
} // namespace ceres

15
extern/libmv/third_party/ceres/internal/ceres/coordinate_descent_minimizer.h vendored
View File

@ -37,12 +37,14 @@
#include "ceres/evaluator.h"
#include "ceres/minimizer.h"
#include "ceres/problem_impl.h"
#include "ceres/program.h"
#include "ceres/solver.h"
namespace ceres {
namespace internal {
class Program;
class LinearSolver;
// Given a Program, and a ParameterBlockOrdering which partitions
// (non-exhaustively) the Hessian matrix into independent sets,
// perform coordinate descent on the parameter blocks in the
@ -66,6 +68,17 @@ class CoordinateDescentMinimizer : public Minimizer {
double* parameters,
Solver::Summary* summary);
// Verify that each group in the ordering forms an independent set.
static bool IsOrderingValid(const Program& program,
const ParameterBlockOrdering& ordering,
string* message);
// Find a recursive decomposition of the Hessian matrix as a set
// of independent sets of decreasing size and invert it. This
// seems to work better in practice, i.e., Cameras before
// points.
static ParameterBlockOrdering* CreateOrdering(const Program& program);
private:
void Solve(Program* program,
LinearSolver* linear_solver,

367
extern/libmv/third_party/ceres/internal/ceres/covariance_impl.cc vendored
View File

@ -38,6 +38,25 @@
#include <cstdlib>
#include <utility>
#include <vector>
#include "Eigen/SparseCore"
// Suppress unused local variable warning from Eigen Ordering.h #included by
// SparseQR in Eigen 3.2.0. This was fixed in Eigen 3.2.1, but 3.2.0 is still
// widely used (Ubuntu 14.04), and Ceres won't compile otherwise due to -Werror.
#if defined(_MSC_VER)
#pragma warning( push )
#pragma warning( disable : 4189 )
#else
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-but-set-variable"
#endif
#include "Eigen/SparseQR"
#if defined(_MSC_VER)
#pragma warning( pop )
#else
#pragma GCC diagnostic pop
#endif
#include "Eigen/SVD"
#include "ceres/compressed_col_sparse_matrix_utils.h"
#include "ceres/compressed_row_sparse_matrix.h"
@ -53,40 +72,6 @@
namespace ceres {
namespace internal {
namespace {
// Per thread storage for SuiteSparse.
#ifndef CERES_NO_SUITESPARSE
struct PerThreadContext {
explicit PerThreadContext(int num_rows)
: solution(NULL),
solution_set(NULL),
y_workspace(NULL),
e_workspace(NULL),
rhs(NULL) {
rhs = ss.CreateDenseVector(NULL, num_rows, num_rows);
}
~PerThreadContext() {
ss.Free(solution);
ss.Free(solution_set);
ss.Free(y_workspace);
ss.Free(e_workspace);
ss.Free(rhs);
}
cholmod_dense* solution;
cholmod_sparse* solution_set;
cholmod_dense* y_workspace;
cholmod_dense* e_workspace;
cholmod_dense* rhs;
SuiteSparse ss;
};
#endif
} // namespace
typedef vector<pair<const double*, const double*> > CovarianceBlocks;
@ -94,8 +79,17 @@ CovarianceImpl::CovarianceImpl(const Covariance::Options& options)
: options_(options),
is_computed_(false),
is_valid_(false) {
evaluate_options_.num_threads = options.num_threads;
evaluate_options_.apply_loss_function = options.apply_loss_function;
#ifndef CERES_USE_OPENMP
if (options_.num_threads > 1) {
LOG(WARNING)
<< "OpenMP support is not compiled into this binary; "
<< "only options.num_threads = 1 is supported. Switching "
<< "to single threaded mode.";
options_.num_threads = 1;
}
#endif
evaluate_options_.num_threads = options_.num_threads;
evaluate_options_.apply_loss_function = options_.apply_loss_function;
}
CovarianceImpl::~CovarianceImpl() {
@ -396,11 +390,15 @@ bool CovarianceImpl::ComputeCovarianceValues() {
case DENSE_SVD:
return ComputeCovarianceValuesUsingDenseSVD();
#ifndef CERES_NO_SUITESPARSE
case SPARSE_CHOLESKY:
return ComputeCovarianceValuesUsingSparseCholesky();
case SPARSE_QR:
return ComputeCovarianceValuesUsingSparseQR();
case SUITE_SPARSE_QR:
return ComputeCovarianceValuesUsingSuiteSparseQR();
#else
LOG(ERROR) << "SuiteSparse is required to use the "
<< "SUITE_SPARSE_QR algorithm.";
return false;
#endif
case EIGEN_SPARSE_QR:
return ComputeCovarianceValuesUsingEigenSparseQR();
default:
LOG(ERROR) << "Unsupported covariance estimation algorithm type: "
<< CovarianceAlgorithmTypeToString(options_.algorithm_type);
@ -409,197 +407,7 @@ bool CovarianceImpl::ComputeCovarianceValues() {
return false;
}
bool CovarianceImpl::ComputeCovarianceValuesUsingSparseCholesky() {
EventLogger event_logger(
"CovarianceImpl::ComputeCovarianceValuesUsingSparseCholesky");
#ifndef CERES_NO_SUITESPARSE
if (covariance_matrix_.get() == NULL) {
// Nothing to do, all zeros covariance matrix.
return true;
}
SuiteSparse ss;
CRSMatrix jacobian;
problem_->Evaluate(evaluate_options_, NULL, NULL, NULL, &jacobian);
event_logger.AddEvent("Evaluate");
// m is a transposed view of the Jacobian.
cholmod_sparse cholmod_jacobian_view;
cholmod_jacobian_view.nrow = jacobian.num_cols;
cholmod_jacobian_view.ncol = jacobian.num_rows;
cholmod_jacobian_view.nzmax = jacobian.values.size();
cholmod_jacobian_view.nz = NULL;
cholmod_jacobian_view.p = reinterpret_cast<void*>(&jacobian.rows[0]);
cholmod_jacobian_view.i = reinterpret_cast<void*>(&jacobian.cols[0]);
cholmod_jacobian_view.x = reinterpret_cast<void*>(&jacobian.values[0]);
cholmod_jacobian_view.z = NULL;
cholmod_jacobian_view.stype = 0; // Matrix is not symmetric.
cholmod_jacobian_view.itype = CHOLMOD_INT;
cholmod_jacobian_view.xtype = CHOLMOD_REAL;
cholmod_jacobian_view.dtype = CHOLMOD_DOUBLE;
cholmod_jacobian_view.sorted = 1;
cholmod_jacobian_view.packed = 1;
string message;
cholmod_factor* factor = ss.AnalyzeCholesky(&cholmod_jacobian_view, &message);
event_logger.AddEvent("Symbolic Factorization");
if (factor == NULL) {
LOG(ERROR) << "Covariance estimation failed. "
<< "CHOLMOD symbolic cholesky factorization returned with: "
<< message;
return false;
}
LinearSolverTerminationType termination_type =
ss.Cholesky(&cholmod_jacobian_view, factor, &message);
event_logger.AddEvent("Numeric Factorization");
if (termination_type != LINEAR_SOLVER_SUCCESS) {
LOG(ERROR) << "Covariance estimation failed. "
<< "CHOLMOD numeric cholesky factorization returned with: "
<< message;
ss.Free(factor);
return false;
}
const double reciprocal_condition_number =
cholmod_rcond(factor, ss.mutable_cc());
if (reciprocal_condition_number <
options_.min_reciprocal_condition_number) {
LOG(ERROR) << "Cholesky factorization of J'J is not reliable. "
<< "Reciprocal condition number: "
<< reciprocal_condition_number << " "
<< "min_reciprocal_condition_number: "
<< options_.min_reciprocal_condition_number;
ss.Free(factor);
return false;
}
const int num_rows = covariance_matrix_->num_rows();
const int* rows = covariance_matrix_->rows();
const int* cols = covariance_matrix_->cols();
double* values = covariance_matrix_->mutable_values();
// The following loop exploits the fact that the i^th column of A^{-1}
// is given by the solution to the linear system
//
// A x = e_i
//
// where e_i is a vector with e(i) = 1 and all other entries zero.
//
// Since the covariance matrix is symmetric, the i^th row and column
// are equal.
//
// The ifdef separates two different version of SuiteSparse. Newer
// versions of SuiteSparse have the cholmod_solve2 function which
// re-uses memory across calls.
#if (SUITESPARSE_VERSION < 4002)
cholmod_dense* rhs = ss.CreateDenseVector(NULL, num_rows, num_rows);
double* rhs_x = reinterpret_cast<double*>(rhs->x);
for (int r = 0; r < num_rows; ++r) {
int row_begin = rows[r];
int row_end = rows[r + 1];
if (row_end == row_begin) {
continue;
}
rhs_x[r] = 1.0;
cholmod_dense* solution = ss.Solve(factor, rhs, &message);
double* solution_x = reinterpret_cast<double*>(solution->x);
for (int idx = row_begin; idx < row_end; ++idx) {
const int c = cols[idx];
values[idx] = solution_x[c];
}
ss.Free(solution);
rhs_x[r] = 0.0;
}
ss.Free(rhs);
#else // SUITESPARSE_VERSION < 4002
const int num_threads = options_.num_threads;
vector<PerThreadContext*> contexts(num_threads);
for (int i = 0; i < num_threads; ++i) {
contexts[i] = new PerThreadContext(num_rows);
}
// The first call to cholmod_solve2 is not thread safe, since it
// changes the factorization from supernodal to simplicial etc.
{
PerThreadContext* context = contexts[0];
double* context_rhs_x = reinterpret_cast<double*>(context->rhs->x);
context_rhs_x[0] = 1.0;
cholmod_solve2(CHOLMOD_A,
factor,
context->rhs,
NULL,
&context->solution,
&context->solution_set,
&context->y_workspace,
&context->e_workspace,
context->ss.mutable_cc());
context_rhs_x[0] = 0.0;
}
#pragma omp parallel for num_threads(num_threads) schedule(dynamic)
for (int r = 0; r < num_rows; ++r) {
int row_begin = rows[r];
int row_end = rows[r + 1];
if (row_end == row_begin) {
continue;
}
# ifdef CERES_USE_OPENMP
int thread_id = omp_get_thread_num();
# else
int thread_id = 0;
# endif
PerThreadContext* context = contexts[thread_id];
double* context_rhs_x = reinterpret_cast<double*>(context->rhs->x);
context_rhs_x[r] = 1.0;
// TODO(sameeragarwal) There should be a more efficient way
// involving the use of Bset but I am unable to make it work right
// now.
cholmod_solve2(CHOLMOD_A,
factor,
context->rhs,
NULL,
&context->solution,
&context->solution_set,
&context->y_workspace,
&context->e_workspace,
context->ss.mutable_cc());
double* solution_x = reinterpret_cast<double*>(context->solution->x);
for (int idx = row_begin; idx < row_end; ++idx) {
const int c = cols[idx];
values[idx] = solution_x[c];
}
context_rhs_x[r] = 0.0;
}
for (int i = 0; i < num_threads; ++i) {
delete contexts[i];
}
#endif // SUITESPARSE_VERSION < 4002
ss.Free(factor);
event_logger.AddEvent("Inversion");
return true;
#else // CERES_NO_SUITESPARSE
return false;
#endif // CERES_NO_SUITESPARSE
};
bool CovarianceImpl::ComputeCovarianceValuesUsingSparseQR() {
bool CovarianceImpl::ComputeCovarianceValuesUsingSuiteSparseQR() {
EventLogger event_logger(
"CovarianceImpl::ComputeCovarianceValuesUsingSparseQR");
@ -851,7 +659,102 @@ bool CovarianceImpl::ComputeCovarianceValuesUsingDenseSVD() {
}
event_logger.AddEvent("CopyToCovarianceMatrix");
return true;
};
}
bool CovarianceImpl::ComputeCovarianceValuesUsingEigenSparseQR() {
EventLogger event_logger(
"CovarianceImpl::ComputeCovarianceValuesUsingEigenSparseQR");
if (covariance_matrix_.get() == NULL) {
// Nothing to do, all zeros covariance matrix.
return true;
}
CRSMatrix jacobian;
problem_->Evaluate(evaluate_options_, NULL, NULL, NULL, &jacobian);
event_logger.AddEvent("Evaluate");
typedef Eigen::SparseMatrix<double, Eigen::ColMajor> EigenSparseMatrix;
// Convert the matrix to column major order as required by SparseQR.
EigenSparseMatrix sparse_jacobian =
Eigen::MappedSparseMatrix<double, Eigen::RowMajor>(
jacobian.num_rows, jacobian.num_cols,
static_cast<int>(jacobian.values.size()),
jacobian.rows.data(), jacobian.cols.data(), jacobian.values.data());
event_logger.AddEvent("ConvertToSparseMatrix");
Eigen::SparseQR<EigenSparseMatrix, Eigen::COLAMDOrdering<int> >
qr_solver(sparse_jacobian);
event_logger.AddEvent("QRDecomposition");
if(qr_solver.info() != Eigen::Success) {
LOG(ERROR) << "Eigen::SparseQR decomposition failed.";
return false;
}
if (qr_solver.rank() < jacobian.num_cols) {
LOG(ERROR) << "Jacobian matrix is rank deficient. "
<< "Number of columns: " << jacobian.num_cols
<< " rank: " << qr_solver.rank();
return false;
}
const int* rows = covariance_matrix_->rows();
const int* cols = covariance_matrix_->cols();
double* values = covariance_matrix_->mutable_values();
// Compute the inverse column permutation used by QR factorization.
Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> inverse_permutation =
qr_solver.colsPermutation().inverse();
// The following loop exploits the fact that the i^th column of A^{-1}
// is given by the solution to the linear system
//
// A x = e_i
//
// where e_i is a vector with e(i) = 1 and all other entries zero.
//
// Since the covariance matrix is symmetric, the i^th row and column
// are equal.
const int num_cols = jacobian.num_cols;
const int num_threads = options_.num_threads;
scoped_array<double> workspace(new double[num_threads * num_cols]);
#pragma omp parallel for num_threads(num_threads) schedule(dynamic)
for (int r = 0; r < num_cols; ++r) {
const int row_begin = rows[r];
const int row_end = rows[r + 1];
if (row_end == row_begin) {
continue;
}
# ifdef CERES_USE_OPENMP
int thread_id = omp_get_thread_num();
# else
int thread_id = 0;
# endif
double* solution = workspace.get() + thread_id * num_cols;
SolveRTRWithSparseRHS<int>(
num_cols,
qr_solver.matrixR().innerIndexPtr(),
qr_solver.matrixR().outerIndexPtr(),
&qr_solver.matrixR().data().value(0),
inverse_permutation.indices().coeff(r),
solution);
// Assign the values of the computed covariance using the
// inverse permutation used in the QR factorization.
for (int idx = row_begin; idx < row_end; ++idx) {
const int c = cols[idx];
values[idx] = solution[inverse_permutation.indices().coeff(c)];
}
}
event_logger.AddEvent("Inverse");
return true;
}
} // namespace internal
} // namespace ceres

4
extern/libmv/third_party/ceres/internal/ceres/covariance_impl.h vendored
View File

@ -64,9 +64,9 @@ class CovarianceImpl {
ProblemImpl* problem);
bool ComputeCovarianceValues();
bool ComputeCovarianceValuesUsingSparseCholesky();
bool ComputeCovarianceValuesUsingSparseQR();
bool ComputeCovarianceValuesUsingDenseSVD();
bool ComputeCovarianceValuesUsingSuiteSparseQR();
bool ComputeCovarianceValuesUsingEigenSparseQR();
const CompressedRowSparseMatrix* covariance_matrix() const {
return covariance_matrix_.get();

7
extern/libmv/third_party/ceres/internal/ceres/cxsparse.h vendored
View File

@ -38,7 +38,6 @@
#include <vector>
#include "cs.h"
#include "ceres/internal/port.h"
namespace ceres {
namespace internal {
@ -130,9 +129,13 @@ class CXSparse {
#else // CERES_NO_CXSPARSE
class CXSparse {};
typedef void cs_dis;
class CXSparse {
public:
void Free(void*) {};
};
#endif // CERES_NO_CXSPARSE
#endif // CERES_INTERNAL_CXSPARSE_H_

11
extern/libmv/third_party/ceres/internal/ceres/dynamic_compressed_row_jacobian_writer.cc vendored
View File

@ -54,8 +54,15 @@ SparseMatrix* DynamicCompressedRowJacobianWriter::CreateJacobian() const {
num_effective_parameters,
0);
CompressedRowJacobianWriter::PopulateJacobianRowAndColumnBlockVectors(
program_, jacobian);
vector<int>* row_blocks = jacobian->mutable_row_blocks();
for (int i = 0; i < jacobian->num_rows(); ++i) {
row_blocks->push_back(1);
}
vector<int>* col_blocks = jacobian->mutable_col_blocks();
for (int i = 0; i < jacobian->num_cols(); ++i) {
col_blocks->push_back(1);
}
return jacobian;
}

11
extern/libmv/third_party/ceres/internal/ceres/gradient_checking_cost_function.cc vendored
View File

@ -310,6 +310,17 @@ ProblemImpl* CreateGradientCheckingProblemImpl(ProblemImpl* problem_impl,
parameter_blocks);
}
// Normally, when a problem is given to the solver, we guarantee
// that the state pointers for each parameter block point to the
// user provided data. Since we are creating this new problem from a
// problem given to us at an arbitrary stage of the solve, we cannot
// depend on this being the case, so we explicitly call
// SetParameterBlockStatePtrsToUserStatePtrs to ensure that this is
// the case.
gradient_checking_problem_impl
->mutable_program()
->SetParameterBlockStatePtrsToUserStatePtrs();
return gradient_checking_problem_impl;
}

81
extern/libmv/third_party/ceres/internal/ceres/gradient_problem.cc vendored Normal file
View File

@ -0,0 +1,81 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#include "ceres/gradient_problem.h"
#include "ceres/local_parameterization.h"
#include "glog/logging.h"
namespace ceres {
GradientProblem::GradientProblem(FirstOrderFunction* function)
: function_(function),
parameterization_(
new IdentityParameterization(function_->NumParameters())),
scratch_(new double[function_->NumParameters()]) {
}
GradientProblem::GradientProblem(FirstOrderFunction* function,
LocalParameterization* parameterization)
: function_(function),
parameterization_(parameterization),
scratch_(new double[function_->NumParameters()]) {
CHECK_EQ(function_->NumParameters(), parameterization_->GlobalSize());
}
int GradientProblem::NumParameters() const {
return function_->NumParameters();
}
int GradientProblem::NumLocalParameters() const {
return parameterization_->LocalSize();
}
bool GradientProblem::Evaluate(const double* parameters,
double* cost,
double* gradient) const {
if (gradient == NULL) {
return function_->Evaluate(parameters, cost, NULL);
}
return (function_->Evaluate(parameters, cost, scratch_.get()) &&
parameterization_->MultiplyByJacobian(parameters,
1,
scratch_.get(),
gradient));
}
bool GradientProblem::Plus(const double* x,
const double* delta,
double* x_plus_delta) const {
return parameterization_->Plus(x, delta, x_plus_delta);
}
} // namespace ceres

98
extern/libmv/third_party/ceres/internal/ceres/gradient_problem_evaluator.h vendored Normal file
View File

@ -0,0 +1,98 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2013 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#ifndef CERES_INTERNAL_GRADIENT_PROBLEM_EVALUATOR_H_
#define CERES_INTERNAL_GRADIENT_PROBLEM_EVALUATOR_H_
#include <map>
#include <string>
#include "ceres/evaluator.h"
#include "ceres/execution_summary.h"
#include "ceres/gradient_problem.h"
#include "ceres/internal/port.h"
#include "ceres/wall_time.h"
namespace ceres {
namespace internal {
class GradientProblemEvaluator : public Evaluator {
public:
explicit GradientProblemEvaluator(const GradientProblem& problem)
: problem_(problem) {}
virtual ~GradientProblemEvaluator() {}
virtual SparseMatrix* CreateJacobian() const { return NULL; }
virtual bool Evaluate(const EvaluateOptions& evaluate_options,
const double* state,
double* cost,
double* residuals,
double* gradient,
SparseMatrix* jacobian) {
CHECK(jacobian == NULL);
ScopedExecutionTimer total_timer("Evaluator::Total", &execution_summary_);
ScopedExecutionTimer call_type_timer(
gradient == NULL ? "Evaluator::Cost" : "Evaluator::Gradient",
&execution_summary_);
return problem_.Evaluate(state, cost, gradient);
}
virtual bool Plus(const double* state,
const double* delta,
double* state_plus_delta) const {
return problem_.Plus(state, delta, state_plus_delta);
}
virtual int NumParameters() const {
return problem_.NumParameters();
}
virtual int NumEffectiveParameters() const {
return problem_.NumLocalParameters();
}
virtual int NumResiduals() const { return 1; }
virtual map<string, int> CallStatistics() const {
return execution_summary_.calls();
}
virtual map<string, double> TimeStatistics() const {
return execution_summary_.times();
}
private:
const GradientProblem& problem_;
::ceres::internal::ExecutionSummary execution_summary_;
};
} // namespace internal
} // namespace ceres
#endif // CERES_INTERNAL_GRADIENT_PROBLEM_EVALUATOR_H_

270
extern/libmv/third_party/ceres/internal/ceres/gradient_problem_solver.cc vendored Normal file
View File

@ -0,0 +1,270 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#include "ceres/gradient_problem_solver.h"
#include "ceres/callbacks.h"
#include "ceres/gradient_problem.h"
#include "ceres/gradient_problem_evaluator.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/port.h"
#include "ceres/map_util.h"
#include "ceres/minimizer.h"
#include "ceres/solver.h"
#include "ceres/solver_utils.h"
#include "ceres/stringprintf.h"
#include "ceres/types.h"
#include "ceres/wall_time.h"
namespace ceres {
using internal::StringPrintf;
using internal::StringAppendF;
namespace {
Solver::Options GradientProblemSolverOptionsToSolverOptions(
const GradientProblemSolver::Options& options) {
#define COPY_OPTION(x) solver_options.x = options.x
Solver::Options solver_options;
solver_options.minimizer_type = LINE_SEARCH;
COPY_OPTION(line_search_direction_type);
COPY_OPTION(line_search_type);
COPY_OPTION(nonlinear_conjugate_gradient_type);
COPY_OPTION(max_lbfgs_rank);
COPY_OPTION(use_approximate_eigenvalue_bfgs_scaling);
COPY_OPTION(line_search_interpolation_type);
COPY_OPTION(min_line_search_step_size);
COPY_OPTION(line_search_sufficient_function_decrease);
COPY_OPTION(max_line_search_step_contraction);
COPY_OPTION(min_line_search_step_contraction);
COPY_OPTION(max_num_line_search_step_size_iterations);
COPY_OPTION(max_num_line_search_direction_restarts);
COPY_OPTION(line_search_sufficient_curvature_decrease);
COPY_OPTION(max_line_search_step_expansion);
COPY_OPTION(max_num_iterations);
COPY_OPTION(max_solver_time_in_seconds);
COPY_OPTION(function_tolerance);
COPY_OPTION(gradient_tolerance);
COPY_OPTION(logging_type);
COPY_OPTION(minimizer_progress_to_stdout);
COPY_OPTION(callbacks);
return solver_options;
#undef COPY_OPTION
}
} // namespace
GradientProblemSolver::~GradientProblemSolver() {
}
void GradientProblemSolver::Solve(const GradientProblemSolver::Options& options,
const GradientProblem& problem,
double* parameters_ptr,
GradientProblemSolver::Summary* summary) {
using internal::scoped_ptr;
using internal::WallTimeInSeconds;
using internal::Minimizer;
using internal::GradientProblemEvaluator;
using internal::LoggingCallback;
using internal::SetSummaryFinalCost;
double start_time = WallTimeInSeconds();
Solver::Options solver_options =
GradientProblemSolverOptionsToSolverOptions(options);
*CHECK_NOTNULL(summary) = Summary();
summary->num_parameters = problem.NumParameters();
summary->num_local_parameters = problem.NumLocalParameters();
summary->line_search_direction_type = options.line_search_direction_type; // NOLINT
summary->line_search_interpolation_type = options.line_search_interpolation_type; // NOLINT
summary->line_search_type = options.line_search_type;
summary->max_lbfgs_rank = options.max_lbfgs_rank;
summary->nonlinear_conjugate_gradient_type = options.nonlinear_conjugate_gradient_type; // NOLINT
// Check validity
if (!solver_options.IsValid(&summary->message)) {
LOG(ERROR) << "Terminating: " << summary->message;
return;
}
// Assuming that the parameter blocks in the program have been
Minimizer::Options minimizer_options;
minimizer_options = Minimizer::Options(solver_options);
minimizer_options.evaluator.reset(new GradientProblemEvaluator(problem));
scoped_ptr<IterationCallback> logging_callback;
if (options.logging_type != SILENT) {
logging_callback.reset(
new LoggingCallback(LINE_SEARCH, options.minimizer_progress_to_stdout));
minimizer_options.callbacks.insert(minimizer_options.callbacks.begin(),
logging_callback.get());
}
scoped_ptr<Minimizer> minimizer(Minimizer::Create(LINE_SEARCH));
Vector solution(problem.NumParameters());
VectorRef parameters(parameters_ptr, problem.NumParameters());
solution = parameters;
Solver::Summary solver_summary;
solver_summary.fixed_cost = 0.0;
solver_summary.preprocessor_time_in_seconds = 0.0;
solver_summary.postprocessor_time_in_seconds = 0.0;
minimizer->Minimize(minimizer_options, solution.data(), &solver_summary);
summary->termination_type = solver_summary.termination_type;
summary->message = solver_summary.message;
summary->initial_cost = solver_summary.initial_cost;
summary->final_cost = solver_summary.final_cost;
summary->iterations = solver_summary.iterations;
if (summary->IsSolutionUsable()) {
parameters = solution;
SetSummaryFinalCost(summary);
}
const map<string, double>& evaluator_time_statistics =
minimizer_options.evaluator->TimeStatistics();
summary->cost_evaluation_time_in_seconds =
FindWithDefault(evaluator_time_statistics, "Evaluator::Residual", 0.0);
summary->gradient_evaluation_time_in_seconds =
FindWithDefault(evaluator_time_statistics, "Evaluator::Jacobian", 0.0);
summary->total_time_in_seconds = WallTimeInSeconds() - start_time;
}
// Invalid values for most fields, to ensure that we are not
// accidentally reporting default values.
GradientProblemSolver::Summary::Summary()
: termination_type(FAILURE),
message("ceres::GradientProblemSolve was not called."),
initial_cost(-1.0),
final_cost(-1.0),
total_time_in_seconds(-1.0),
cost_evaluation_time_in_seconds(-1.0),
gradient_evaluation_time_in_seconds(-1.0),
num_parameters(-1),
num_local_parameters(-1),
line_search_direction_type(LBFGS),
line_search_type(ARMIJO),
line_search_interpolation_type(BISECTION),
nonlinear_conjugate_gradient_type(FLETCHER_REEVES),
max_lbfgs_rank(-1) {
}
bool GradientProblemSolver::Summary::IsSolutionUsable() const {
return internal::IsSolutionUsable(*this);
}
string GradientProblemSolver::Summary::BriefReport() const {
return StringPrintf("Ceres GradientProblemSolver Report: "
"Iterations: %d, "
"Initial cost: %e, "
"Final cost: %e, "
"Termination: %s",
static_cast<int>(iterations.size()),
initial_cost,
final_cost,
TerminationTypeToString(termination_type));
};
string GradientProblemSolver::Summary::FullReport() const {
using internal::VersionString;
string report = string("\nSolver Summary (v " + VersionString() + ")\n\n");
StringAppendF(&report, "Parameters % 25d\n", num_parameters);
if (num_local_parameters != num_parameters) {
StringAppendF(&report, "Local parameters % 25d\n",
num_local_parameters);
}
string line_search_direction_string;
if (line_search_direction_type == LBFGS) {
line_search_direction_string = StringPrintf("LBFGS (%d)", max_lbfgs_rank);
} else if (line_search_direction_type == NONLINEAR_CONJUGATE_GRADIENT) {
line_search_direction_string =
NonlinearConjugateGradientTypeToString(
nonlinear_conjugate_gradient_type);
} else {
line_search_direction_string =
LineSearchDirectionTypeToString(line_search_direction_type);
}
StringAppendF(&report, "Line search direction %19s\n",
line_search_direction_string.c_str());
const string line_search_type_string =
StringPrintf("%s %s",
LineSearchInterpolationTypeToString(
line_search_interpolation_type),
LineSearchTypeToString(line_search_type));
StringAppendF(&report, "Line search type %19s\n",
line_search_type_string.c_str());
StringAppendF(&report, "\n");
StringAppendF(&report, "\nCost:\n");
StringAppendF(&report, "Initial % 30e\n", initial_cost);
if (termination_type != FAILURE &&
termination_type != USER_FAILURE) {
StringAppendF(&report, "Final % 30e\n", final_cost);
StringAppendF(&report, "Change % 30e\n",
initial_cost - final_cost);
}
StringAppendF(&report, "\nMinimizer iterations % 16d\n",
static_cast<int>(iterations.size()));
StringAppendF(&report, "\nTime (in seconds):\n");
StringAppendF(&report, "\n Cost evaluation %23.3f\n",
cost_evaluation_time_in_seconds);
StringAppendF(&report, " Gradient evaluation %23.3f\n",
gradient_evaluation_time_in_seconds);
StringAppendF(&report, "Total %25.3f\n\n",
total_time_in_seconds);
StringAppendF(&report, "Termination: %25s (%s)\n",
TerminationTypeToString(termination_type), message.c_str());
return report;
};
void Solve(const GradientProblemSolver::Options& options,
const GradientProblem& problem,
double* parameters,
GradientProblemSolver::Summary* summary) {
GradientProblemSolver solver;
solver.Solve(options, problem, parameters, summary);
}
} // namespace ceres

68
extern/libmv/third_party/ceres/internal/ceres/graph.h vendored
View File

@ -43,13 +43,75 @@
namespace ceres {
namespace internal {
// A weighted undirected graph templated over the vertex ids. Vertex
// should be hashable and comparable.
// A unweighted undirected graph templated over the vertex ids. Vertex
// should be hashable.
template <typename Vertex>
class Graph {
public:
Graph() {}
// Add a vertex.
void AddVertex(const Vertex& vertex) {
if (vertices_.insert(vertex).second) {
edges_[vertex] = HashSet<Vertex>();
}
}
bool RemoveVertex(const Vertex& vertex) {
if (vertices_.find(vertex) == vertices_.end()) {
return false;
}
vertices_.erase(vertex);
const HashSet<Vertex>& sinks = edges_[vertex];
for (typename HashSet<Vertex>::const_iterator it = sinks.begin();
it != sinks.end(); ++it) {
edges_[*it].erase(vertex);
}
edges_.erase(vertex);
return true;
}
// Add an edge between the vertex1 and vertex2. Calling AddEdge on a
// pair of vertices which do not exist in the graph yet will result
// in undefined behavior.
//
// It is legal to call this method repeatedly for the same set of
// vertices.
void AddEdge(const Vertex& vertex1, const Vertex& vertex2) {
DCHECK(vertices_.find(vertex1) != vertices_.end());
DCHECK(vertices_.find(vertex2) != vertices_.end());
if (edges_[vertex1].insert(vertex2).second) {
edges_[vertex2].insert(vertex1);
}
}
// Calling Neighbors on a vertex not in the graph will result in
// undefined behaviour.
const HashSet<Vertex>& Neighbors(const Vertex& vertex) const {
return FindOrDie(edges_, vertex);
}
const HashSet<Vertex>& vertices() const {
return vertices_;
}
private:
HashSet<Vertex> vertices_;
HashMap<Vertex, HashSet<Vertex> > edges_;
CERES_DISALLOW_COPY_AND_ASSIGN(Graph);
};
// A weighted undirected graph templated over the vertex ids. Vertex
// should be hashable and comparable.
template <typename Vertex>
class WeightedGraph {
public:
WeightedGraph() {}
// Add a weighted vertex. If the vertex already exists in the graph,
// its weight is set to the new weight.
void AddVertex(const Vertex& vertex, double weight) {
@ -152,7 +214,7 @@ class Graph {
HashMap<Vertex, HashSet<Vertex> > edges_;
HashMap<pair<Vertex, Vertex>, double> edge_weights_;
CERES_DISALLOW_COPY_AND_ASSIGN(Graph);
CERES_DISALLOW_COPY_AND_ASSIGN(WeightedGraph);
};
} // namespace internal

7
extern/libmv/third_party/ceres/internal/ceres/graph_algorithms.h vendored
View File

@ -38,6 +38,7 @@
#include <utility>
#include "ceres/collections_port.h"
#include "ceres/graph.h"
#include "ceres/wall_time.h"
#include "glog/logging.h"
namespace ceres {
@ -270,11 +271,11 @@ Vertex FindConnectedComponent(const Vertex& vertex,
// spanning forest, or a collection of linear paths that span the
// graph G.
template <typename Vertex>
Graph<Vertex>*
Degree2MaximumSpanningForest(const Graph<Vertex>& graph) {
WeightedGraph<Vertex>*
Degree2MaximumSpanningForest(const WeightedGraph<Vertex>& graph) {
// Array of edges sorted in decreasing order of their weights.
vector<pair<double, pair<Vertex, Vertex> > > weighted_edges;
Graph<Vertex>* forest = new Graph<Vertex>();
WeightedGraph<Vertex>* forest = new WeightedGraph<Vertex>();
// Disjoint-set to keep track of the connected components in the
// maximum spanning tree.

1
extern/libmv/third_party/ceres/internal/ceres/iterative_schur_complement_solver.cc vendored
View File

@ -101,6 +101,7 @@ LinearSolver::Summary IterativeSchurComplementSolver::SolveImpl(
// complement matrix with the block diagonal of the matrix F'F as
// the preconditioner.
LinearSolver::Options cg_options;
cg_options.min_num_iterations = options_.min_num_iterations;
cg_options.max_num_iterations = options_.max_num_iterations;
ConjugateGradientsSolver cg_solver(cg_options);
LinearSolver::PerSolveOptions cg_per_solve_options;

2
extern/libmv/third_party/ceres/internal/ceres/line_search_direction.cc vendored
View File

@ -65,7 +65,7 @@ class NonlinearConjugateGradient : public LineSearchDirection {
case FLETCHER_REEVES:
beta = current.gradient_squared_norm / previous.gradient_squared_norm;
break;
case POLAK_RIBIRERE:
case POLAK_RIBIERE:
gradient_change = current.gradient - previous.gradient;
beta = (current.gradient.dot(gradient_change) /
previous.gradient_squared_norm);

5
extern/libmv/third_party/ceres/internal/ceres/line_search_minimizer.cc vendored
View File

@ -103,7 +103,7 @@ void LineSearchMinimizer::Minimize(const Minimizer::Options& options,
double start_time = WallTimeInSeconds();
double iteration_start_time = start_time;
Evaluator* evaluator = CHECK_NOTNULL(options.evaluator);
Evaluator* evaluator = CHECK_NOTNULL(options.evaluator.get());
const int num_parameters = evaluator->NumParameters();
const int num_effective_parameters = evaluator->NumEffectiveParameters();
@ -375,7 +375,6 @@ void LineSearchMinimizer::Minimize(const Minimizer::Options& options,
WallTimeInSeconds() - start_time
+ summary->preprocessor_time_in_seconds;
summary->iterations.push_back(iteration_summary);
++summary->num_successful_steps;
if (iteration_summary.gradient_max_norm <= options.gradient_tolerance) {
@ -401,6 +400,8 @@ void LineSearchMinimizer::Minimize(const Minimizer::Options& options,
VLOG_IF(1, is_not_silent) << "Terminating: " << summary->message;
break;
}
summary->iterations.push_back(iteration_summary);
}
}

106
extern/libmv/third_party/ceres/internal/ceres/line_search_preprocessor.cc vendored Normal file
View File

@ -0,0 +1,106 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#include "ceres/line_search_preprocessor.h"
#include <numeric>
#include <string>
#include "ceres/evaluator.h"
#include "ceres/minimizer.h"
#include "ceres/problem_impl.h"
#include "ceres/program.h"
#include "ceres/wall_time.h"
namespace ceres {
namespace internal {
namespace {
bool IsProgramValid(const Program& program, string* error) {
if (program.IsBoundsConstrained()) {
*error = "LINE_SEARCH Minimizer does not support bounds.";
return false;
}
return program.ParameterBlocksAreFinite(error);
}
bool SetupEvaluator(PreprocessedProblem* pp) {
pp->evaluator_options = Evaluator::Options();
// This ensures that we get a Block Jacobian Evaluator without any
// requirement on orderings.
pp->evaluator_options.linear_solver_type = CGNR;
pp->evaluator_options.num_eliminate_blocks = 0;
pp->evaluator_options.num_threads = pp->options.num_threads;
pp->evaluator.reset(Evaluator::Create(pp->evaluator_options,
pp->reduced_program.get(),
&pp->error));
return (pp->evaluator.get() != NULL);
}
} // namespace
LineSearchPreprocessor::~LineSearchPreprocessor() {
}
bool LineSearchPreprocessor::Preprocess(const Solver::Options& options,
ProblemImpl* problem,
PreprocessedProblem* pp) {
CHECK_NOTNULL(pp);
pp->options = options;
ChangeNumThreadsIfNeeded(&pp->options);
pp->problem = problem;
Program* program = problem->mutable_program();
if (!IsProgramValid(*program, &pp->error)) {
return false;
}
pp->reduced_program.reset(
program->CreateReducedProgram(&pp->removed_parameter_blocks,
&pp->fixed_cost,
&pp->error));
if (pp->reduced_program.get() == NULL) {
return false;
}
if (pp->reduced_program->NumParameterBlocks() == 0) {
return true;
}
if (!SetupEvaluator(pp)) {
return false;
}
SetupCommonMinimizerOptions(pp);
return true;
}
} // namespace internal
} // namespace ceres

50
extern/libmv/third_party/ceres/internal/ceres/line_search_preprocessor.h vendored Normal file
View File

@ -0,0 +1,50 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 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: sameragarwal@google.com (Sameer Agarwal)
#ifndef CERES_INTERNAL_LINE_SEARCH_PREPROCESSOR_H_
#define CERES_INTERNAL_LINE_SEARCH_PREPROCESSOR_H_
#include "ceres/preprocessor.h"
namespace ceres {
namespace internal {
class LineSearchPreprocessor : public Preprocessor {
public:
virtual ~LineSearchPreprocessor();
virtual bool Preprocess(const Solver::Options& options,
ProblemImpl* problem,
PreprocessedProblem* preprocessed_problem);
};
} // namespace internal
} // namespace ceres
#endif // CERES_INTERNAL_LINE_SEARCH_PREPROCESSOR_H_

38
extern/libmv/third_party/ceres/internal/ceres/linear_solver.cc vendored
View File

@ -45,26 +45,48 @@ namespace internal {
LinearSolver::~LinearSolver() {
}
LinearSolverType LinearSolver::LinearSolverForZeroEBlocks(
LinearSolverType linear_solver_type) {
if (!IsSchurType(linear_solver_type)) {
return linear_solver_type;
}
if (linear_solver_type == SPARSE_SCHUR) {
return SPARSE_NORMAL_CHOLESKY;
}
if (linear_solver_type == DENSE_SCHUR) {
// TODO(sameeragarwal): This is probably not a great choice.
// Ideally, we should have a DENSE_NORMAL_CHOLESKY, that can take
// a BlockSparseMatrix as input.
return DENSE_QR;
}
if (linear_solver_type == ITERATIVE_SCHUR) {
return CGNR;
}
return linear_solver_type;
}
LinearSolver* LinearSolver::Create(const LinearSolver::Options& options) {
switch (options.type) {
case CGNR:
return new CgnrSolver(options);
case SPARSE_NORMAL_CHOLESKY:
#if defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CXSPARSE)
LOG(WARNING) << "SPARSE_NORMAL_CHOLESKY is not available. Please "
<< "build Ceres with SuiteSparse or CXSparse. "
<< "Returning NULL.";
#if defined(CERES_NO_SUITESPARSE) && \
defined(CERES_NO_CXSPARSE) && \
!defined(CERES_USE_EIGEN_SPARSE)
return NULL;
#else
return new SparseNormalCholeskySolver(options);
#endif
case SPARSE_SCHUR:
#if defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CXSPARSE)
LOG(WARNING) << "SPARSE_SCHUR is not available. Please "
<< "build Ceres with SuiteSparse or CXSparse. "
<< "Returning NULL.";
#if defined(CERES_NO_SUITESPARSE) && \
defined(CERES_NO_CXSPARSE) && \
!defined(CERES_USE_EIGEN_SPARSE)
return NULL;
#else
return new SparseSchurComplementSolver(options);

8
extern/libmv/third_party/ceres/internal/ceres/linear_solver.h vendored
View File

@ -274,6 +274,14 @@ class LinearSolver {
string message;
};
// If the optimization problem is such that there are no remaining
// e-blocks, a Schur type linear solver cannot be used. If the
// linear solver is of Schur type, this function implements a policy
// to select an alternate nearest linear solver to the one selected
// by the user. The input linear_solver_type is returned otherwise.
static LinearSolverType LinearSolverForZeroEBlocks(
LinearSolverType linear_solver_type);
virtual ~LinearSolver();
// Solve Ax = b.

44
extern/libmv/third_party/ceres/internal/ceres/local_parameterization.cc vendored
View File

@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
@ -36,6 +36,23 @@
namespace ceres {
LocalParameterization::~LocalParameterization() {
}
bool LocalParameterization::MultiplyByJacobian(const double* x,
const int num_rows,
const double* global_matrix,
double* local_matrix) const {
Matrix jacobian(GlobalSize(), LocalSize());
if (!ComputeJacobian(x, jacobian.data())) {
return false;
}
MatrixRef(local_matrix, num_rows, LocalSize()) =
ConstMatrixRef(global_matrix, num_rows, GlobalSize()) * jacobian;
return true;
}
IdentityParameterization::IdentityParameterization(const int size)
: size_(size) {
CHECK_GT(size, 0);
@ -55,6 +72,16 @@ bool IdentityParameterization::ComputeJacobian(const double* x,
return true;
}
bool IdentityParameterization::MultiplyByJacobian(const double* x,
const int num_cols,
const double* global_matrix,
double* local_matrix) const {
std::copy(global_matrix,
global_matrix + num_cols * GlobalSize(),
local_matrix);
return true;
}
SubsetParameterization::SubsetParameterization(
int size,
const vector<int>& constant_parameters)
@ -108,6 +135,21 @@ bool SubsetParameterization::ComputeJacobian(const double* x,
return true;
}
bool SubsetParameterization::MultiplyByJacobian(const double* x,
const int num_rows,
const double* global_matrix,
double* local_matrix) const {
for (int row = 0; row < num_rows; ++row) {
for (int col = 0, j = 0; col < constancy_mask_.size(); ++col) {
if (!constancy_mask_[col]) {
local_matrix[row * LocalSize() + j++] =
global_matrix[row * GlobalSize() + col];
}
}
}
return true;
}
bool QuaternionParameterization::Plus(const double* x,
const double* delta,
double* x_plus_delta) const {

35
extern/libmv/third_party/ceres/internal/ceres/loss_function.cc vendored
View File

@ -34,13 +34,14 @@
#include <cmath>
#include <cstddef>
#include <limits>
namespace ceres {
void TrivialLoss::Evaluate(double s, double rho[3]) const {
rho[0] = s;
rho[1] = 1;
rho[2] = 0;
rho[1] = 1.0;
rho[2] = 0.0;
}
void HuberLoss::Evaluate(double s, double rho[3]) const {
@ -48,32 +49,32 @@ void HuberLoss::Evaluate(double s, double rho[3]) const {
// Outlier region.
// 'r' is always positive.
const double r = sqrt(s);
rho[0] = 2 * a_ * r - b_;
rho[1] = a_ / r;
rho[2] = - rho[1] / (2 * s);
rho[0] = 2.0 * a_ * r - b_;
rho[1] = std::max(std::numeric_limits<double>::min(), a_ / r);
rho[2] = - rho[1] / (2.0 * s);
} else {
// Inlier region.
rho[0] = s;
rho[1] = 1;
rho[2] = 0;
rho[1] = 1.0;
rho[2] = 0.0;
}
}
void SoftLOneLoss::Evaluate(double s, double rho[3]) const {
const double sum = 1 + s * c_;
const double sum = 1.0 + s * c_;
const double tmp = sqrt(sum);
// 'sum' and 'tmp' are always positive, assuming that 's' is.
rho[0] = 2 * b_ * (tmp - 1);
rho[1] = 1 / tmp;
rho[2] = - (c_ * rho[1]) / (2 * sum);
rho[0] = 2.0 * b_ * (tmp - 1.0);
rho[1] = std::max(std::numeric_limits<double>::min(), 1.0 / tmp);
rho[2] = - (c_ * rho[1]) / (2.0 * sum);
}
void CauchyLoss::Evaluate(double s, double rho[3]) const {
const double sum = 1 + s * c_;
const double inv = 1 / sum;
const double sum = 1.0 + s * c_;
const double inv = 1.0 / sum;
// 'sum' and 'inv' are always positive, assuming that 's' is.
rho[0] = b_ * log(sum);
rho[1] = inv;
rho[1] = std::max(std::numeric_limits<double>::min(), inv);
rho[2] = - c_ * (inv * inv);
}
@ -82,8 +83,8 @@ void ArctanLoss::Evaluate(double s, double rho[3]) const {
const double inv = 1 / sum;
// 'sum' and 'inv' are always positive.
rho[0] = a_ * atan2(s, a_);
rho[1] = inv;
rho[2] = -2 * s * b_ * (inv * inv);
rho[1] = std::max(std::numeric_limits<double>::min(), inv);
rho[2] = -2.0 * s * b_ * (inv * inv);
}
TolerantLoss::TolerantLoss(double a, double b)
@ -108,7 +109,7 @@ void TolerantLoss::Evaluate(double s, double rho[3]) const {
} else {
const double e_x = exp(x);
rho[0] = b_ * log(1.0 + e_x) - c_;
rho[1] = e_x / (1.0 + e_x);
rho[1] = std::max(std::numeric_limits<double>::min(), e_x / (1.0 + e_x));
rho[2] = 0.5 / (b_ * (1.0 + cosh(x)));
}
}

16
extern/libmv/third_party/ceres/internal/ceres/minimizer.cc vendored
View File

@ -28,13 +28,29 @@
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#include "ceres/line_search_minimizer.h"
#include "ceres/minimizer.h"
#include "ceres/trust_region_minimizer.h"
#include "ceres/types.h"
#include "glog/logging.h"
namespace ceres {
namespace internal {
Minimizer* Minimizer::Create(MinimizerType minimizer_type) {
if (minimizer_type == TRUST_REGION) {
return new TrustRegionMinimizer;
}
if (minimizer_type == LINE_SEARCH) {
return new LineSearchMinimizer;
}
LOG(FATAL) << "Unknown minimizer_type: " << minimizer_type;
return NULL;
}
Minimizer::~Minimizer() {}
bool Minimizer::RunCallbacks(const Minimizer::Options& options,

35
extern/libmv/third_party/ceres/internal/ceres/minimizer.h vendored
View File

@ -41,9 +41,10 @@ namespace ceres {
namespace internal {
class Evaluator;
class LinearSolver;
class SparseMatrix;
class TrustRegionStrategy;
class CoordinateDescentMinimizer;
class LinearSolver;
// Interface for non-linear least squares solvers.
class Minimizer {
@ -107,14 +108,10 @@ class Minimizer {
options.line_search_sufficient_curvature_decrease;
max_line_search_step_expansion =
options.max_line_search_step_expansion;
is_silent = (options.logging_type == SILENT);
evaluator = NULL;
trust_region_strategy = NULL;
jacobian = NULL;
callbacks = options.callbacks;
inner_iteration_minimizer = NULL;
inner_iteration_tolerance = options.inner_iteration_tolerance;
is_silent = (options.logging_type == SILENT);
is_constrained = false;
callbacks = options.callbacks;
}
int max_num_iterations;
@ -154,10 +151,14 @@ class Minimizer {
int max_num_line_search_direction_restarts;
double line_search_sufficient_curvature_decrease;
double max_line_search_step_expansion;
double inner_iteration_tolerance;
// If true, then all logging is disabled.
bool is_silent;
// Use a bounds constrained optimization algorithm.
bool is_constrained;
// List of callbacks that are executed by the Minimizer at the end
// of each iteration.
//
@ -165,27 +166,23 @@ class Minimizer {
vector<IterationCallback*> callbacks;
// Object responsible for evaluating the cost, residuals and
// Jacobian matrix. The Options struct does not own this pointer.
Evaluator* evaluator;
// Jacobian matrix.
shared_ptr<Evaluator> evaluator;
// Object responsible for actually computing the trust region
// step, and sizing the trust region radius. The Options struct
// does not own this pointer.
TrustRegionStrategy* trust_region_strategy;
// step, and sizing the trust region radius.
shared_ptr<TrustRegionStrategy> trust_region_strategy;
// Object holding the Jacobian matrix. It is assumed that the
// sparsity structure of the matrix has already been initialized
// and will remain constant for the life time of the
// optimization. The Options struct does not own this pointer.
SparseMatrix* jacobian;
// optimization.
shared_ptr<SparseMatrix> jacobian;
Minimizer* inner_iteration_minimizer;
double inner_iteration_tolerance;
// Use a bounds constrained optimization algorithm.
bool is_constrained;
shared_ptr<CoordinateDescentMinimizer> inner_iteration_minimizer;
};
static Minimizer* Create(MinimizerType minimizer_type);
static bool RunCallbacks(const Options& options,
const IterationSummary& iteration_summary,
Solver::Summary* summary);

27
extern/libmv/third_party/ceres/internal/ceres/parameter_block_ordering.cc vendored
View File

@ -37,6 +37,7 @@
#include "ceres/parameter_block.h"
#include "ceres/program.h"
#include "ceres/residual_block.h"
#include "ceres/wall_time.h"
#include "glog/logging.h"
namespace ceres {
@ -45,8 +46,10 @@ namespace internal {
int ComputeStableSchurOrdering(const Program& program,
vector<ParameterBlock*>* ordering) {
CHECK_NOTNULL(ordering)->clear();
EventLogger event_logger("ComputeStableSchurOrdering");
scoped_ptr<Graph< ParameterBlock*> > graph(CreateHessianGraph(program));
event_logger.AddEvent("CreateHessianGraph");
const vector<ParameterBlock*>& parameter_blocks = program.parameter_blocks();
const HashSet<ParameterBlock*>& vertices = graph->vertices();
for (int i = 0; i < parameter_blocks.size(); ++i) {
@ -54,8 +57,10 @@ int ComputeStableSchurOrdering(const Program& program,
ordering->push_back(parameter_blocks[i]);
}
}
event_logger.AddEvent("Preordering");
int independent_set_size = StableIndependentSetOrdering(*graph, ordering);
event_logger.AddEvent("StableIndependentSet");
// Add the excluded blocks to back of the ordering vector.
for (int i = 0; i < parameter_blocks.size(); ++i) {
@ -64,6 +69,7 @@ int ComputeStableSchurOrdering(const Program& program,
ordering->push_back(parameter_block);
}
}
event_logger.AddEvent("ConstantParameterBlocks");
return independent_set_size;
}
@ -109,8 +115,7 @@ void ComputeRecursiveIndependentSetOrdering(const Program& program,
}
}
Graph<ParameterBlock*>*
CreateHessianGraph(const Program& program) {
Graph<ParameterBlock*>* CreateHessianGraph(const Program& program) {
Graph<ParameterBlock*>* graph = CHECK_NOTNULL(new Graph<ParameterBlock*>);
const vector<ParameterBlock*>& parameter_blocks = program.parameter_blocks();
for (int i = 0; i < parameter_blocks.size(); ++i) {
@ -144,5 +149,21 @@ CreateHessianGraph(const Program& program) {
return graph;
}
void OrderingToGroupSizes(const ParameterBlockOrdering* ordering,
vector<int>* group_sizes) {
CHECK_NOTNULL(group_sizes)->clear();
if (ordering == NULL) {
return;
}
const map<int, set<double*> >& group_to_elements =
ordering->group_to_elements();
for (map<int, set<double*> >::const_iterator it = group_to_elements.begin();
it != group_to_elements.end();
++it) {
group_sizes->push_back(it->second.size());
}
}
} // namespace internal
} // namespace ceres

5
extern/libmv/third_party/ceres/internal/ceres/parameter_block_ordering.h vendored
View File

@ -78,6 +78,11 @@ void ComputeRecursiveIndependentSetOrdering(const Program& program,
// parameter blocks, if they co-occur in a residual block.
Graph<ParameterBlock*>* CreateHessianGraph(const Program& program);
// Iterate over each of the groups in order of their priority and fill
// summary with their sizes.
void OrderingToGroupSizes(const ParameterBlockOrdering* ordering,
vector<int>* group_sizes);
} // namespace internal
} // namespace ceres

10
extern/libmv/third_party/ceres/internal/ceres/preconditioner.cc vendored
View File

@ -37,6 +37,16 @@ namespace internal {
Preconditioner::~Preconditioner() {
}
PreconditionerType Preconditioner::PreconditionerForZeroEBlocks(
PreconditionerType preconditioner_type) {
if (preconditioner_type == SCHUR_JACOBI ||
preconditioner_type == CLUSTER_JACOBI ||
preconditioner_type == CLUSTER_TRIDIAGONAL) {
return JACOBI;
}
return preconditioner_type;
}
SparseMatrixPreconditionerWrapper::SparseMatrixPreconditionerWrapper(
const SparseMatrix* matrix)
: matrix_(CHECK_NOTNULL(matrix)) {

9
extern/libmv/third_party/ceres/internal/ceres/preconditioner.h vendored
View File

@ -36,6 +36,7 @@
#include "ceres/compressed_row_sparse_matrix.h"
#include "ceres/linear_operator.h"
#include "ceres/sparse_matrix.h"
#include "ceres/types.h"
namespace ceres {
namespace internal {
@ -95,6 +96,14 @@ class Preconditioner : public LinearOperator {
int f_block_size;
};
// If the optimization problem is such that there are no remaining
// e-blocks, ITERATIVE_SCHUR with a Schur type preconditioner cannot
// be used. This function returns JACOBI if a preconditioner for
// ITERATIVE_SCHUR is used. The input preconditioner_type is
// returned otherwise.
static PreconditionerType PreconditionerForZeroEBlocks(
PreconditionerType preconditioner_type);
virtual ~Preconditioner();
// Update the numerical value of the preconditioner for the linear

113
extern/libmv/third_party/ceres/internal/ceres/preprocessor.cc vendored Normal file
View File

@ -0,0 +1,113 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 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: sameragarwal@google.com (Sameer Agarwal)
#include "ceres/callbacks.h"
#include "ceres/gradient_checking_cost_function.h"
#include "ceres/line_search_preprocessor.h"
#include "ceres/preprocessor.h"
#include "ceres/problem_impl.h"
#include "ceres/solver.h"
#include "ceres/trust_region_preprocessor.h"
namespace ceres {
namespace internal {
Preprocessor* Preprocessor::Create(MinimizerType minimizer_type) {
if (minimizer_type == TRUST_REGION) {
return new TrustRegionPreprocessor;
}
if (minimizer_type == LINE_SEARCH) {
return new LineSearchPreprocessor;
}
LOG(FATAL) << "Unknown minimizer_type: " << minimizer_type;
return NULL;
}
Preprocessor::~Preprocessor() {
}
void ChangeNumThreadsIfNeeded(Solver::Options* options) {
#ifndef CERES_USE_OPENMP
if (options->num_threads > 1) {
LOG(WARNING)
<< "OpenMP support is not compiled into this binary; "
<< "only options.num_threads = 1 is supported. Switching "
<< "to single threaded mode.";
options->num_threads = 1;
}
// Only the Trust Region solver currently uses a linear solver.
if (options->minimizer_type == TRUST_REGION &&
options->num_linear_solver_threads > 1) {
LOG(WARNING)
<< "OpenMP support is not compiled into this binary; "
<< "only options.num_linear_solver_threads=1 is supported. Switching "
<< "to single threaded mode.";
options->num_linear_solver_threads = 1;
}
#endif // CERES_USE_OPENMP
}
void SetupCommonMinimizerOptions(PreprocessedProblem* pp) {
const Solver::Options& options = pp->options;
Program* program = pp->reduced_program.get();
// Assuming that the parameter blocks in the program have been
// reordered as needed, extract them into a contiguous vector.
pp->reduced_parameters.resize(program->NumParameters());
double* reduced_parameters = pp->reduced_parameters.data();
program->ParameterBlocksToStateVector(reduced_parameters);
Minimizer::Options& minimizer_options = pp->minimizer_options;
minimizer_options = Minimizer::Options(options);
minimizer_options.evaluator = pp->evaluator;
if (options.logging_type != SILENT) {
pp->logging_callback.reset(
new LoggingCallback(options.minimizer_type,
options.minimizer_progress_to_stdout));
minimizer_options.callbacks.insert(minimizer_options.callbacks.begin(),
pp->logging_callback.get());
}
if (options.update_state_every_iteration) {
pp->state_updating_callback.reset(
new StateUpdatingCallback(program, reduced_parameters));
// This must get pushed to the front of the callbacks so that it
// is run before any of the user callbacks.
minimizer_options.callbacks.insert(minimizer_options.callbacks.begin(),
pp->state_updating_callback.get());
}
}
} // namespace internal
} // namespace ceres

119
extern/libmv/third_party/ceres/internal/ceres/preprocessor.h vendored Normal file
View File

@ -0,0 +1,119 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 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: sameragarwal@google.com (Sameer Agarwal)
#ifndef CERES_INTERNAL_PREPROCESSOR_H_
#define CERES_INTERNAL_PREPROCESSOR_H_
#include "ceres/coordinate_descent_minimizer.h"
#include "ceres/evaluator.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/port.h"
#include "ceres/internal/scoped_ptr.h"
#include "ceres/iteration_callback.h"
#include "ceres/linear_solver.h"
#include "ceres/minimizer.h"
#include "ceres/problem_impl.h"
#include "ceres/program.h"
#include "ceres/solver.h"
namespace ceres {
namespace internal {
struct PreprocessedProblem;
// Given a Problem object and a Solver::Options object indicating the
// configuration of the solver, the job of the Preprocessor is to
// analyze the Problem and perform the setup needed to solve it using
// the desired Minimization algorithm. The setup involves removing
// redundancies in the input problem (inactive parameter and residual
// blocks), finding fill reducing orderings as needed, configuring and
// creating various objects needed by the Minimizer to solve the
// problem such as an evaluator, a linear solver etc.
//
// Each Minimizer (LineSearchMinimizer and TrustRegionMinimizer) comes
// with a corresponding Preprocessor (LineSearchPreprocessor and
// TrustRegionPreprocessor) that knows about its needs and performs
// the preprocessing needed.
//
// The output of the Preprocessor is stored in a PreprocessedProblem
// object.
class Preprocessor {
public:
// Factory.
static Preprocessor* Create(MinimizerType minimizer_type);
virtual ~Preprocessor();
virtual bool Preprocess(const Solver::Options& options,
ProblemImpl* problem,
PreprocessedProblem* pp) = 0;
};
// A PreprocessedProblem is the result of running the Preprocessor on
// a Problem and Solver::Options object.
struct PreprocessedProblem {
PreprocessedProblem()
: fixed_cost(0.0) {
}
string error;
Solver::Options options;
LinearSolver::Options linear_solver_options;
Evaluator::Options evaluator_options;
Minimizer::Options minimizer_options;
ProblemImpl* problem;
scoped_ptr<ProblemImpl> gradient_checking_problem;
scoped_ptr<Program> reduced_program;
scoped_ptr<LinearSolver> linear_solver;
scoped_ptr<IterationCallback> logging_callback;
scoped_ptr<IterationCallback> state_updating_callback;
shared_ptr<Evaluator> evaluator;
shared_ptr<CoordinateDescentMinimizer> inner_iteration_minimizer;
vector<double*> removed_parameter_blocks;
Vector reduced_parameters;
double fixed_cost;
};
// Common functions used by various preprocessors.
// If OpenMP support is not available and user has requested more than
// one thread, then set the *_num_threads options as needed to 1.
void ChangeNumThreadsIfNeeded(Solver::Options* options);
// Extract the effective parameter vector from the preprocessed
// problem and setup bits of the Minimizer::Options object that are
// common to all Preprocessors.
void SetupCommonMinimizerOptions(PreprocessedProblem* pp);
} // namespace internal
} // namespace ceres
#endif // CERES_INTERNAL_PREPROCESSOR_H_

10
extern/libmv/third_party/ceres/internal/ceres/problem.cc vendored
View File

@ -251,6 +251,16 @@ void Problem::GetParameterBlocksForResidualBlock(
parameter_blocks);
}
const CostFunction* Problem::GetCostFunctionForResidualBlock(
const ResidualBlockId residual_block) const {
return problem_impl_->GetCostFunctionForResidualBlock(residual_block);
}
const LossFunction* Problem::GetLossFunctionForResidualBlock(
const ResidualBlockId residual_block) const {
return problem_impl_->GetLossFunctionForResidualBlock(residual_block);
}
void Problem::GetResidualBlocksForParameterBlock(
const double* values,
vector<ResidualBlockId>* residual_blocks) const {

10
extern/libmv/third_party/ceres/internal/ceres/problem_impl.cc vendored
View File

@ -823,6 +823,16 @@ void ProblemImpl::GetParameterBlocksForResidualBlock(
}
}
const CostFunction* ProblemImpl::GetCostFunctionForResidualBlock(
const ResidualBlockId residual_block) const {
return residual_block->cost_function();
}
const LossFunction* ProblemImpl::GetLossFunctionForResidualBlock(
const ResidualBlockId residual_block) const {
return residual_block->loss_function();
}
void ProblemImpl::GetResidualBlocksForParameterBlock(
const double* values,
vector<ResidualBlockId>* residual_blocks) const {

5
extern/libmv/third_party/ceres/internal/ceres/problem_impl.h vendored
View File

@ -157,6 +157,11 @@ class ProblemImpl {
const ResidualBlockId residual_block,
vector<double*>* parameter_blocks) const;
const CostFunction* GetCostFunctionForResidualBlock(
const ResidualBlockId residual_block) const;
const LossFunction* GetLossFunctionForResidualBlock(
const ResidualBlockId residual_block) const;
void GetResidualBlocksForParameterBlock(
const double* values,
vector<ResidualBlockId>* residual_blocks) const;

255
extern/libmv/third_party/ceres/internal/ceres/program.cc vendored
View File

@ -32,6 +32,7 @@
#include <map>
#include <vector>
#include "ceres/array_utils.h"
#include "ceres/casts.h"
#include "ceres/compressed_row_sparse_matrix.h"
#include "ceres/cost_function.h"
@ -44,6 +45,7 @@
#include "ceres/problem.h"
#include "ceres/residual_block.h"
#include "ceres/stl_util.h"
#include "ceres/triplet_sparse_matrix.h"
namespace ceres {
namespace internal {
@ -170,6 +172,259 @@ bool Program::IsValid() const {
return true;
}
bool Program::ParameterBlocksAreFinite(string* message) const {
CHECK_NOTNULL(message);
for (int i = 0; i < parameter_blocks_.size(); ++i) {
const ParameterBlock* parameter_block = parameter_blocks_[i];
const double* array = parameter_block->user_state();
const int size = parameter_block->Size();
const int invalid_index = FindInvalidValue(size, array);
if (invalid_index != size) {
*message = StringPrintf(
"ParameterBlock: %p with size %d has at least one invalid value.\n"
"First invalid value is at index: %d.\n"
"Parameter block values: ",
array, size, invalid_index);
AppendArrayToString(size, array, message);
return false;
}
}
return true;
}
bool Program::IsBoundsConstrained() const {
for (int i = 0; i < parameter_blocks_.size(); ++i) {
const ParameterBlock* parameter_block = parameter_blocks_[i];
if (parameter_block->IsConstant()) {
continue;
}
const int size = parameter_block->Size();
for (int j = 0; j < size; ++j) {
const double lower_bound = parameter_block->LowerBoundForParameter(j);
const double upper_bound = parameter_block->UpperBoundForParameter(j);
if (lower_bound > -std::numeric_limits<double>::max() ||
upper_bound < std::numeric_limits<double>::max()) {
return true;
}
}
}
return false;
}
bool Program::IsFeasible(string* message) const {
CHECK_NOTNULL(message);
for (int i = 0; i < parameter_blocks_.size(); ++i) {
const ParameterBlock* parameter_block = parameter_blocks_[i];
const double* parameters = parameter_block->user_state();
const int size = parameter_block->Size();
if (parameter_block->IsConstant()) {
// Constant parameter blocks must start in the feasible region
// to ultimately produce a feasible solution, since Ceres cannot
// change them.
for (int j = 0; j < size; ++j) {
const double lower_bound = parameter_block->LowerBoundForParameter(j);
const double upper_bound = parameter_block->UpperBoundForParameter(j);
if (parameters[j] < lower_bound || parameters[j] > upper_bound) {
*message = StringPrintf(
"ParameterBlock: %p with size %d has at least one infeasible "
"value."
"\nFirst infeasible value is at index: %d."
"\nLower bound: %e, value: %e, upper bound: %e"
"\nParameter block values: ",
parameters, size, j, lower_bound, parameters[j], upper_bound);
AppendArrayToString(size, parameters, message);
return false;
}
}
} else {
// Variable parameter blocks must have non-empty feasible
// regions, otherwise there is no way to produce a feasible
// solution.
for (int j = 0; j < size; ++j) {
const double lower_bound = parameter_block->LowerBoundForParameter(j);
const double upper_bound = parameter_block->UpperBoundForParameter(j);
if (lower_bound >= upper_bound) {
*message = StringPrintf(
"ParameterBlock: %p with size %d has at least one infeasible "
"bound."
"\nFirst infeasible bound is at index: %d."
"\nLower bound: %e, upper bound: %e"
"\nParameter block values: ",
parameters, size, j, lower_bound, upper_bound);
AppendArrayToString(size, parameters, message);
return false;
}
}
}
}
return true;
}
Program* Program::CreateReducedProgram(vector<double*>* removed_parameter_blocks,
double* fixed_cost,
string* error) const {
CHECK_NOTNULL(removed_parameter_blocks);
CHECK_NOTNULL(fixed_cost);
CHECK_NOTNULL(error);
scoped_ptr<Program> reduced_program(new Program(*this));
if (!reduced_program->RemoveFixedBlocks(removed_parameter_blocks,
fixed_cost,
error)) {
return NULL;
}
reduced_program->SetParameterOffsetsAndIndex();
return reduced_program.release();
}
bool Program::RemoveFixedBlocks(vector<double*>* removed_parameter_blocks,
double* fixed_cost,
string* error) {
CHECK_NOTNULL(removed_parameter_blocks);
CHECK_NOTNULL(fixed_cost);
CHECK_NOTNULL(error);
scoped_array<double> residual_block_evaluate_scratch;
residual_block_evaluate_scratch.reset(
new double[MaxScratchDoublesNeededForEvaluate()]);
*fixed_cost = 0.0;
// Mark all the parameters as unused. Abuse the index member of the
// parameter blocks for the marking.
for (int i = 0; i < parameter_blocks_.size(); ++i) {
parameter_blocks_[i]->set_index(-1);
}
// Filter out residual that have all-constant parameters, and mark
// all the parameter blocks that appear in residuals.
int num_active_residual_blocks = 0;
for (int i = 0; i < residual_blocks_.size(); ++i) {
ResidualBlock* residual_block = residual_blocks_[i];
int num_parameter_blocks = residual_block->NumParameterBlocks();
// Determine if the residual block is fixed, and also mark varying
// parameters that appear in the residual block.
bool all_constant = true;
for (int k = 0; k < num_parameter_blocks; k++) {
ParameterBlock* parameter_block = residual_block->parameter_blocks()[k];
if (!parameter_block->IsConstant()) {
all_constant = false;
parameter_block->set_index(1);
}
}
if (!all_constant) {
residual_blocks_[num_active_residual_blocks++] = residual_block;
continue;
}
// The residual is constant and will be removed, so its cost is
// added to the variable fixed_cost.
double cost = 0.0;
if (!residual_block->Evaluate(true,
&cost,
NULL,
NULL,
residual_block_evaluate_scratch.get())) {
*error = StringPrintf("Evaluation of the residual %d failed during "
"removal of fixed residual blocks.", i);
return false;
}
*fixed_cost += cost;
}
residual_blocks_.resize(num_active_residual_blocks);
// Filter out unused or fixed parameter blocks.
int num_active_parameter_blocks = 0;
removed_parameter_blocks->clear();
for (int i = 0; i < parameter_blocks_.size(); ++i) {
ParameterBlock* parameter_block = parameter_blocks_[i];
if (parameter_block->index() == -1) {
removed_parameter_blocks->push_back(parameter_block->mutable_user_state());
} else {
parameter_blocks_[num_active_parameter_blocks++] = parameter_block;
}
}
parameter_blocks_.resize(num_active_parameter_blocks);
if (!(((NumResidualBlocks() == 0) &&
(NumParameterBlocks() == 0)) ||
((NumResidualBlocks() != 0) &&
(NumParameterBlocks() != 0)))) {
*error = "Congratulations, you found a bug in Ceres. Please report it.";
return false;
}
return true;
}
bool Program::IsParameterBlockSetIndependent(const set<double*>& independent_set) const {
// Loop over each residual block and ensure that no two parameter
// blocks in the same residual block are part of
// parameter_block_ptrs as that would violate the assumption that it
// is an independent set in the Hessian matrix.
for (vector<ResidualBlock*>::const_iterator it = residual_blocks_.begin();
it != residual_blocks_.end();
++it) {
ParameterBlock* const* parameter_blocks = (*it)->parameter_blocks();
const int num_parameter_blocks = (*it)->NumParameterBlocks();
int count = 0;
for (int i = 0; i < num_parameter_blocks; ++i) {
count += independent_set.count(
parameter_blocks[i]->mutable_user_state());
}
if (count > 1) {
return false;
}
}
return true;
}
TripletSparseMatrix* Program::CreateJacobianBlockSparsityTranspose() const {
// Matrix to store the block sparsity structure of the Jacobian.
TripletSparseMatrix* tsm =
new TripletSparseMatrix(NumParameterBlocks(),
NumResidualBlocks(),
10 * NumResidualBlocks());
int num_nonzeros = 0;
int* rows = tsm->mutable_rows();
int* cols = tsm->mutable_cols();
double* values = tsm->mutable_values();
for (int c = 0; c < residual_blocks_.size(); ++c) {
const ResidualBlock* residual_block = residual_blocks_[c];
const int num_parameter_blocks = residual_block->NumParameterBlocks();
ParameterBlock* const* parameter_blocks =
residual_block->parameter_blocks();
for (int j = 0; j < num_parameter_blocks; ++j) {
if (parameter_blocks[j]->IsConstant()) {
continue;
}
// Re-size the matrix if needed.
if (num_nonzeros >= tsm->max_num_nonzeros()) {
tsm->set_num_nonzeros(num_nonzeros);
tsm->Reserve(2 * num_nonzeros);
rows = tsm->mutable_rows();
cols = tsm->mutable_cols();
values = tsm->mutable_values();
}
const int r = parameter_blocks[j]->index();
rows[num_nonzeros] = r;
cols[num_nonzeros] = c;
values[num_nonzeros] = 1.0;
++num_nonzeros;
}
}
tsm->set_num_nonzeros(num_nonzeros);
return tsm;
}
int Program::NumResidualBlocks() const {
return residual_blocks_.size();
}

58
extern/libmv/third_party/ceres/internal/ceres/program.h vendored
View File

@ -31,6 +31,7 @@
#ifndef CERES_INTERNAL_PROGRAM_H_
#define CERES_INTERNAL_PROGRAM_H_
#include <set>
#include <string>
#include <vector>
#include "ceres/internal/port.h"
@ -41,6 +42,7 @@ namespace internal {
class ParameterBlock;
class ProblemImpl;
class ResidualBlock;
class TripletSparseMatrix;
// A nonlinear least squares optimization problem. This is different from the
// similarly-named "Problem" object, which offers a mutation interface for
@ -103,6 +105,47 @@ class Program {
// offsets) are correct.
bool IsValid() const;
bool ParameterBlocksAreFinite(string* message) const;
// Returns true if the program has any non-constant parameter blocks
// which have non-trivial bounds constraints.
bool IsBoundsConstrained() const;
// Returns false, if the program has any constant parameter blocks
// which are not feasible, or any variable parameter blocks which
// have a lower bound greater than or equal to the upper bound.
bool IsFeasible(string* message) const;
// Loop over each residual block and ensure that no two parameter
// blocks in the same residual block are part of
// parameter_blocks as that would violate the assumption that it
// is an independent set in the Hessian matrix.
bool IsParameterBlockSetIndependent(const set<double*>& independent_set) const;
// Create a TripletSparseMatrix which contains the zero-one
// structure corresponding to the block sparsity of the transpose of
// the Jacobian matrix.
//
// Caller owns the result.
TripletSparseMatrix* CreateJacobianBlockSparsityTranspose() const;
// Create a copy of this program and removes constant parameter
// blocks and residual blocks with no varying parameter blocks while
// preserving their relative order.
//
// removed_parameter_blocks on exit will contain the list of
// parameter blocks that were removed.
//
// fixed_cost will be equal to the sum of the costs of the residual
// blocks that were removed.
//
// If there was a problem, then the function will return a NULL
// pointer and error will contain a human readable description of
// the problem.
Program* CreateReducedProgram(vector<double*>* removed_parameter_blocks,
double* fixed_cost,
string* error) const;
// See problem.h for what these do.
int NumParameterBlocks() const;
int NumParameters() const;
@ -120,6 +163,21 @@ class Program {
string ToString() const;
private:
// Remove constant parameter blocks and residual blocks with no
// varying parameter blocks while preserving their relative order.
//
// removed_parameter_blocks on exit will contain the list of
// parameter blocks that were removed.
//
// fixed_cost will be equal to the sum of the costs of the residual
// blocks that were removed.
//
// If there was a problem, then the function will return false and
// error will contain a human readable description of the problem.
bool RemoveFixedBlocks(vector<double*>* removed_parameter_blocks,
double* fixed_cost,
string* message);
// The Program does not own the ParameterBlock or ResidualBlock objects.
vector<ParameterBlock*> parameter_blocks_;
vector<ResidualBlock*> residual_blocks_;

578
extern/libmv/third_party/ceres/internal/ceres/reorder_program.cc vendored Normal file
View File

@ -0,0 +1,578 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#include "ceres/reorder_program.h"
#include <algorithm>
#include <numeric>
#include <vector>
#include "ceres/cxsparse.h"
#include "ceres/internal/port.h"
#include "ceres/ordered_groups.h"
#include "ceres/parameter_block.h"
#include "ceres/parameter_block_ordering.h"
#include "ceres/problem_impl.h"
#include "ceres/program.h"
#include "ceres/residual_block.h"
#include "ceres/solver.h"
#include "ceres/suitesparse.h"
#include "ceres/triplet_sparse_matrix.h"
#include "ceres/types.h"
#include "Eigen/SparseCore"
#ifdef CERES_USE_EIGEN_SPARSE
#include "Eigen/OrderingMethods"
#endif
#include "glog/logging.h"
namespace ceres {
namespace internal {
namespace {
// Find the minimum index of any parameter block to the given
// residual. Parameter blocks that have indices greater than
// size_of_first_elimination_group are considered to have an index
// equal to size_of_first_elimination_group.
static int MinParameterBlock(const ResidualBlock* residual_block,
int size_of_first_elimination_group) {
int min_parameter_block_position = size_of_first_elimination_group;
for (int i = 0; i < residual_block->NumParameterBlocks(); ++i) {
ParameterBlock* parameter_block = residual_block->parameter_blocks()[i];
if (!parameter_block->IsConstant()) {
CHECK_NE(parameter_block->index(), -1)
<< "Did you forget to call Program::SetParameterOffsetsAndIndex()? "
<< "This is a Ceres bug; please contact the developers!";
min_parameter_block_position = std::min(parameter_block->index(),
min_parameter_block_position);
}
}
return min_parameter_block_position;
}
#if EIGEN_VERSION_AT_LEAST(3, 2, 2) && defined(CERES_USE_EIGEN_SPARSE)
Eigen::SparseMatrix<int> CreateBlockJacobian(
const TripletSparseMatrix& block_jacobian_transpose) {
typedef Eigen::SparseMatrix<int> SparseMatrix;
typedef Eigen::Triplet<int> Triplet;
const int* rows = block_jacobian_transpose.rows();
const int* cols = block_jacobian_transpose.cols();
int num_nonzeros = block_jacobian_transpose.num_nonzeros();
std::vector<Triplet> triplets;
triplets.reserve(num_nonzeros);
for (int i = 0; i < num_nonzeros; ++i) {
triplets.push_back(Triplet(cols[i], rows[i], 1));
}
SparseMatrix block_jacobian(block_jacobian_transpose.num_cols(),
block_jacobian_transpose.num_rows());
block_jacobian.setFromTriplets(triplets.begin(), triplets.end());
return block_jacobian;
}
#endif
void OrderingForSparseNormalCholeskyUsingSuiteSparse(
const TripletSparseMatrix& tsm_block_jacobian_transpose,
const vector<ParameterBlock*>& parameter_blocks,
const ParameterBlockOrdering& parameter_block_ordering,
int* ordering) {
#ifdef CERES_NO_SUITESPARSE
LOG(FATAL) << "Congratulations, you found a Ceres bug! "
<< "Please report this error to the developers.";
#else
SuiteSparse ss;
cholmod_sparse* block_jacobian_transpose =
ss.CreateSparseMatrix(
const_cast<TripletSparseMatrix*>(&tsm_block_jacobian_transpose));
// No CAMD or the user did not supply a useful ordering, then just
// use regular AMD.
if (parameter_block_ordering.NumGroups() <= 1 ||
!SuiteSparse::IsConstrainedApproximateMinimumDegreeOrderingAvailable()) {
ss.ApproximateMinimumDegreeOrdering(block_jacobian_transpose, &ordering[0]);
} else {
vector<int> constraints;
for (int i = 0; i < parameter_blocks.size(); ++i) {
constraints.push_back(
parameter_block_ordering.GroupId(
parameter_blocks[i]->mutable_user_state()));
}
ss.ConstrainedApproximateMinimumDegreeOrdering(block_jacobian_transpose,
&constraints[0],
ordering);
}
ss.Free(block_jacobian_transpose);
#endif // CERES_NO_SUITESPARSE
}
void OrderingForSparseNormalCholeskyUsingCXSparse(
const TripletSparseMatrix& tsm_block_jacobian_transpose,
int* ordering) {
#ifdef CERES_NO_CXSPARSE
LOG(FATAL) << "Congratulations, you found a Ceres bug! "
<< "Please report this error to the developers.";
#else // CERES_NO_CXSPARSE
// CXSparse works with J'J instead of J'. So compute the block
// sparsity for J'J and compute an approximate minimum degree
// ordering.
CXSparse cxsparse;
cs_di* block_jacobian_transpose;
block_jacobian_transpose =
cxsparse.CreateSparseMatrix(
const_cast<TripletSparseMatrix*>(&tsm_block_jacobian_transpose));
cs_di* block_jacobian = cxsparse.TransposeMatrix(block_jacobian_transpose);
cs_di* block_hessian =
cxsparse.MatrixMatrixMultiply(block_jacobian_transpose, block_jacobian);
cxsparse.Free(block_jacobian);
cxsparse.Free(block_jacobian_transpose);
cxsparse.ApproximateMinimumDegreeOrdering(block_hessian, ordering);
cxsparse.Free(block_hessian);
#endif // CERES_NO_CXSPARSE
}
#if EIGEN_VERSION_AT_LEAST(3, 2, 2)
void OrderingForSparseNormalCholeskyUsingEigenSparse(
const TripletSparseMatrix& tsm_block_jacobian_transpose,
int* ordering) {
#ifndef CERES_USE_EIGEN_SPARSE
LOG(FATAL) <<
"SPARSE_NORMAL_CHOLESKY cannot be used with EIGEN_SPARSE "
"because Ceres was not built with support for "
"Eigen's SimplicialLDLT decomposition. "
"This requires enabling building with -DEIGENSPARSE=ON.";
#else
// This conversion from a TripletSparseMatrix to a Eigen::Triplet
// matrix is unfortunate, but unavoidable for now. It is not a
// significant performance penalty in the grand scheme of
// things. The right thing to do here would be to get a compressed
// row sparse matrix representation of the jacobian and go from
// there. But that is a project for another day.
typedef Eigen::SparseMatrix<int> SparseMatrix;
const SparseMatrix block_jacobian =
CreateBlockJacobian(tsm_block_jacobian_transpose);
const SparseMatrix block_hessian =
block_jacobian.transpose() * block_jacobian;
Eigen::AMDOrdering<int> amd_ordering;
Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic, int> perm;
amd_ordering(block_hessian, perm);
for (int i = 0; i < block_hessian.rows(); ++i) {
ordering[i] = perm.indices()[i];
}
#endif // CERES_USE_EIGEN_SPARSE
}
#endif
} // namespace
bool ApplyOrdering(const ProblemImpl::ParameterMap& parameter_map,
const ParameterBlockOrdering& ordering,
Program* program,
string* error) {
const int num_parameter_blocks = program->NumParameterBlocks();
if (ordering.NumElements() != num_parameter_blocks) {
*error = StringPrintf("User specified ordering does not have the same "
"number of parameters as the problem. The problem"
"has %d blocks while the ordering has %d blocks.",
num_parameter_blocks,
ordering.NumElements());
return false;
}
vector<ParameterBlock*>* parameter_blocks =
program->mutable_parameter_blocks();
parameter_blocks->clear();
const map<int, set<double*> >& groups =
ordering.group_to_elements();
for (map<int, set<double*> >::const_iterator group_it = groups.begin();
group_it != groups.end();
++group_it) {
const set<double*>& group = group_it->second;
for (set<double*>::const_iterator parameter_block_ptr_it = group.begin();
parameter_block_ptr_it != group.end();
++parameter_block_ptr_it) {
ProblemImpl::ParameterMap::const_iterator parameter_block_it =
parameter_map.find(*parameter_block_ptr_it);
if (parameter_block_it == parameter_map.end()) {
*error = StringPrintf("User specified ordering contains a pointer "
"to a double that is not a parameter block in "
"the problem. The invalid double is in group: %d",
group_it->first);
return false;
}
parameter_blocks->push_back(parameter_block_it->second);
}
}
return true;
}
bool LexicographicallyOrderResidualBlocks(
const int size_of_first_elimination_group,
Program* program,
string* error) {
CHECK_GE(size_of_first_elimination_group, 1)
<< "Congratulations, you found a Ceres bug! Please report this error "
<< "to the developers.";
// Create a histogram of the number of residuals for each E block. There is an
// extra bucket at the end to catch all non-eliminated F blocks.
vector<int> residual_blocks_per_e_block(size_of_first_elimination_group + 1);
vector<ResidualBlock*>* residual_blocks = program->mutable_residual_blocks();
vector<int> min_position_per_residual(residual_blocks->size());
for (int i = 0; i < residual_blocks->size(); ++i) {
ResidualBlock* residual_block = (*residual_blocks)[i];
int position = MinParameterBlock(residual_block,
size_of_first_elimination_group);
min_position_per_residual[i] = position;
DCHECK_LE(position, size_of_first_elimination_group);
residual_blocks_per_e_block[position]++;
}
// Run a cumulative sum on the histogram, to obtain offsets to the start of
// each histogram bucket (where each bucket is for the residuals for that
// E-block).
vector<int> offsets(size_of_first_elimination_group + 1);
std::partial_sum(residual_blocks_per_e_block.begin(),
residual_blocks_per_e_block.end(),
offsets.begin());
CHECK_EQ(offsets.back(), residual_blocks->size())
<< "Congratulations, you found a Ceres bug! Please report this error "
<< "to the developers.";
CHECK(find(residual_blocks_per_e_block.begin(),
residual_blocks_per_e_block.end() - 1, 0) !=
residual_blocks_per_e_block.end())
<< "Congratulations, you found a Ceres bug! Please report this error "
<< "to the developers.";
// Fill in each bucket with the residual blocks for its corresponding E block.
// Each bucket is individually filled from the back of the bucket to the front
// of the bucket. The filling order among the buckets is dictated by the
// residual blocks. This loop uses the offsets as counters; subtracting one
// from each offset as a residual block is placed in the bucket. When the
// filling is finished, the offset pointerts should have shifted down one
// entry (this is verified below).
vector<ResidualBlock*> reordered_residual_blocks(
(*residual_blocks).size(), static_cast<ResidualBlock*>(NULL));
for (int i = 0; i < residual_blocks->size(); ++i) {
int bucket = min_position_per_residual[i];
// Decrement the cursor, which should now point at the next empty position.
offsets[bucket]--;
// Sanity.
CHECK(reordered_residual_blocks[offsets[bucket]] == NULL)
<< "Congratulations, you found a Ceres bug! Please report this error "
<< "to the developers.";
reordered_residual_blocks[offsets[bucket]] = (*residual_blocks)[i];
}
// Sanity check #1: The difference in bucket offsets should match the
// histogram sizes.
for (int i = 0; i < size_of_first_elimination_group; ++i) {
CHECK_EQ(residual_blocks_per_e_block[i], offsets[i + 1] - offsets[i])
<< "Congratulations, you found a Ceres bug! Please report this error "
<< "to the developers.";
}
// Sanity check #2: No NULL's left behind.
for (int i = 0; i < reordered_residual_blocks.size(); ++i) {
CHECK(reordered_residual_blocks[i] != NULL)
<< "Congratulations, you found a Ceres bug! Please report this error "
<< "to the developers.";
}
// Now that the residuals are collected by E block, swap them in place.
swap(*program->mutable_residual_blocks(), reordered_residual_blocks);
return true;
}
// Pre-order the columns corresponding to the schur complement if
// possible.
void MaybeReorderSchurComplementColumnsUsingSuiteSparse(
const ParameterBlockOrdering& parameter_block_ordering,
Program* program) {
#ifndef CERES_NO_SUITESPARSE
SuiteSparse ss;
if (!SuiteSparse::IsConstrainedApproximateMinimumDegreeOrderingAvailable()) {
return;
}
vector<int> constraints;
vector<ParameterBlock*>& parameter_blocks =
*(program->mutable_parameter_blocks());
for (int i = 0; i < parameter_blocks.size(); ++i) {
constraints.push_back(
parameter_block_ordering.GroupId(
parameter_blocks[i]->mutable_user_state()));
}
// Renumber the entries of constraints to be contiguous integers
// as camd requires that the group ids be in the range [0,
// parameter_blocks.size() - 1].
MapValuesToContiguousRange(constraints.size(), &constraints[0]);
// Compute a block sparse presentation of J'.
scoped_ptr<TripletSparseMatrix> tsm_block_jacobian_transpose(
program->CreateJacobianBlockSparsityTranspose());
cholmod_sparse* block_jacobian_transpose =
ss.CreateSparseMatrix(tsm_block_jacobian_transpose.get());
vector<int> ordering(parameter_blocks.size(), 0);
ss.ConstrainedApproximateMinimumDegreeOrdering(block_jacobian_transpose,
&constraints[0],
&ordering[0]);
ss.Free(block_jacobian_transpose);
const vector<ParameterBlock*> parameter_blocks_copy(parameter_blocks);
for (int i = 0; i < program->NumParameterBlocks(); ++i) {
parameter_blocks[i] = parameter_blocks_copy[ordering[i]];
}
program->SetParameterOffsetsAndIndex();
#endif
}
void MaybeReorderSchurComplementColumnsUsingEigen(
const int size_of_first_elimination_group,
const ProblemImpl::ParameterMap& parameter_map,
Program* program) {
#if !EIGEN_VERSION_AT_LEAST(3, 2, 2) || !defined(CERES_USE_EIGEN_SPARSE)
return;
#else
scoped_ptr<TripletSparseMatrix> tsm_block_jacobian_transpose(
program->CreateJacobianBlockSparsityTranspose());
typedef Eigen::SparseMatrix<int> SparseMatrix;
const SparseMatrix block_jacobian =
CreateBlockJacobian(*tsm_block_jacobian_transpose);
const int num_rows = block_jacobian.rows();
const int num_cols = block_jacobian.cols();
// Vertically partition the jacobian in parameter blocks of type E
// and F.
const SparseMatrix E =
block_jacobian.block(0,
0,
num_rows,
size_of_first_elimination_group);
const SparseMatrix F =
block_jacobian.block(0,
size_of_first_elimination_group,
num_rows,
num_cols - size_of_first_elimination_group);
// Block sparsity pattern of the schur complement.
const SparseMatrix block_schur_complement =
F.transpose() * F - F.transpose() * E * E.transpose() * F;
Eigen::AMDOrdering<int> amd_ordering;
Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic, int> perm;
amd_ordering(block_schur_complement, perm);
const vector<ParameterBlock*>& parameter_blocks = program->parameter_blocks();
vector<ParameterBlock*> ordering(num_cols);
// The ordering of the first size_of_first_elimination_group does
// not matter, so we preserve the existing ordering.
for (int i = 0; i < size_of_first_elimination_group; ++i) {
ordering[i] = parameter_blocks[i];
}
// For the rest of the blocks, use the ordering computed using AMD.
for (int i = 0; i < block_schur_complement.cols(); ++i) {
ordering[size_of_first_elimination_group + i] =
parameter_blocks[size_of_first_elimination_group + perm.indices()[i]];
}
swap(*program->mutable_parameter_blocks(), ordering);
program->SetParameterOffsetsAndIndex();
#endif
}
bool ReorderProgramForSchurTypeLinearSolver(
const LinearSolverType linear_solver_type,
const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
const ProblemImpl::ParameterMap& parameter_map,
ParameterBlockOrdering* parameter_block_ordering,
Program* program,
string* error) {
if (parameter_block_ordering->NumElements() !=
program->NumParameterBlocks()) {
*error = StringPrintf(
"The program has %d parameter blocks, but the parameter block "
"ordering has %d parameter blocks.",
program->NumParameterBlocks(),
parameter_block_ordering->NumElements());
return false;
}
if (parameter_block_ordering->NumGroups() == 1) {
// If the user supplied an parameter_block_ordering with just one
// group, it is equivalent to the user supplying NULL as an
// parameter_block_ordering. Ceres is completely free to choose the
// parameter block ordering as it sees fit. For Schur type solvers,
// this means that the user wishes for Ceres to identify the
// e_blocks, which we do by computing a maximal independent set.
vector<ParameterBlock*> schur_ordering;
const int size_of_first_elimination_group =
ComputeStableSchurOrdering(*program, &schur_ordering);
CHECK_EQ(schur_ordering.size(), program->NumParameterBlocks())
<< "Congratulations, you found a Ceres bug! Please report this error "
<< "to the developers.";
// Update the parameter_block_ordering object.
for (int i = 0; i < schur_ordering.size(); ++i) {
double* parameter_block = schur_ordering[i]->mutable_user_state();
const int group_id = (i < size_of_first_elimination_group) ? 0 : 1;
parameter_block_ordering->AddElementToGroup(parameter_block, group_id);
}
// We could call ApplyOrdering but this is cheaper and
// simpler.
swap(*program->mutable_parameter_blocks(), schur_ordering);
} else {
// The user provided an ordering with more than one elimination
// group.
// Verify that the first elimination group is an independent set.
const set<double*>& first_elimination_group =
parameter_block_ordering
->group_to_elements()
.begin()
->second;
if (!program->IsParameterBlockSetIndependent(first_elimination_group)) {
*error =
StringPrintf("The first elimination group in the parameter block "
"ordering of size %zd is not an independent set",
first_elimination_group.size());
return false;
}
if (!ApplyOrdering(parameter_map,
*parameter_block_ordering,
program,
error)) {
return false;
}
}
program->SetParameterOffsetsAndIndex();
const int size_of_first_elimination_group =
parameter_block_ordering->group_to_elements().begin()->second.size();
if (linear_solver_type == SPARSE_SCHUR) {
if (sparse_linear_algebra_library_type == SUITE_SPARSE) {
MaybeReorderSchurComplementColumnsUsingSuiteSparse(
*parameter_block_ordering,
program);
} else if (sparse_linear_algebra_library_type == EIGEN_SPARSE) {
MaybeReorderSchurComplementColumnsUsingEigen(
size_of_first_elimination_group,
parameter_map,
program);
}
}
// Schur type solvers also require that their residual blocks be
// lexicographically ordered.
if (!LexicographicallyOrderResidualBlocks(size_of_first_elimination_group,
program,
error)) {
return false;
}
return true;
}
bool ReorderProgramForSparseNormalCholesky(
const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
const ParameterBlockOrdering& parameter_block_ordering,
Program* program,
string* error) {
// Compute a block sparse presentation of J'.
scoped_ptr<TripletSparseMatrix> tsm_block_jacobian_transpose(
program->CreateJacobianBlockSparsityTranspose());
vector<int> ordering(program->NumParameterBlocks(), 0);
vector<ParameterBlock*>& parameter_blocks =
*(program->mutable_parameter_blocks());
if (sparse_linear_algebra_library_type == SUITE_SPARSE) {
OrderingForSparseNormalCholeskyUsingSuiteSparse(
*tsm_block_jacobian_transpose,
parameter_blocks,
parameter_block_ordering,
&ordering[0]);
} else if (sparse_linear_algebra_library_type == CX_SPARSE) {
OrderingForSparseNormalCholeskyUsingCXSparse(
*tsm_block_jacobian_transpose,
&ordering[0]);
} else if (sparse_linear_algebra_library_type == EIGEN_SPARSE) {
#if EIGEN_VERSION_AT_LEAST(3, 2, 2)
OrderingForSparseNormalCholeskyUsingEigenSparse(
*tsm_block_jacobian_transpose,
&ordering[0]);
#else
// For Eigen versions less than 3.2.2, there is nothing to do as
// older versions of Eigen do not expose a method for doing
// symbolic analysis on pre-ordered matrices, so a block
// pre-ordering is a bit pointless.
return true;
#endif
}
// Apply ordering.
const vector<ParameterBlock*> parameter_blocks_copy(parameter_blocks);
for (int i = 0; i < program->NumParameterBlocks(); ++i) {
parameter_blocks[i] = parameter_blocks_copy[ordering[i]];
}
program->SetParameterOffsetsAndIndex();
return true;
}
} // namespace internal
} // namespace ceres

101
extern/libmv/third_party/ceres/internal/ceres/reorder_program.h vendored Normal file
View File

@ -0,0 +1,101 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#ifndef CERES_INTERNAL_REORDER_PROGRAM_H_
#define CERES_INTERNAL_REORDER_PROGRAM_H_
#include <string>
#include "ceres/internal/port.h"
#include "ceres/parameter_block_ordering.h"
#include "ceres/problem_impl.h"
#include "ceres/types.h"
namespace ceres {
namespace internal {
class Program;
// Reorder the parameter blocks in program using the ordering
bool ApplyOrdering(const ProblemImpl::ParameterMap& parameter_map,
const ParameterBlockOrdering& ordering,
Program* program,
string* error);
// Reorder the residuals for program, if necessary, so that the residuals
// involving each E block occur together. This is a necessary condition for the
// Schur eliminator, which works on these "row blocks" in the jacobian.
bool LexicographicallyOrderResidualBlocks(int size_of_first_elimination_group,
Program* program,
string* error);
// Schur type solvers require that all parameter blocks eliminated
// by the Schur eliminator occur before others and the residuals be
// sorted in lexicographic order of their parameter blocks.
//
// If the parameter_block_ordering only contains one elimination
// group then a maximal independent set is computed and used as the
// first elimination group, otherwise the user's ordering is used.
//
// If the linear solver type is SPARSE_SCHUR and support for
// constrained fill-reducing ordering is available in the sparse
// linear algebra library (SuiteSparse version >= 4.2.0) then
// columns of the schur complement matrix are ordered to reduce the
// fill-in the Cholesky factorization.
//
// Upon return, ordering contains the parameter block ordering that
// was used to order the program.
bool ReorderProgramForSchurTypeLinearSolver(
LinearSolverType linear_solver_type,
SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
const ProblemImpl::ParameterMap& parameter_map,
ParameterBlockOrdering* parameter_block_ordering,
Program* program,
string* error);
// Sparse cholesky factorization routines when doing the sparse
// cholesky factorization of the Jacobian matrix, reorders its
// columns to reduce the fill-in. Compute this permutation and
// re-order the parameter blocks.
//
// When using SuiteSparse, if the parameter_block_ordering contains
// more than one elimination group and support for constrained
// fill-reducing ordering is available in the sparse linear algebra
// library (SuiteSparse version >= 4.2.0) then the fill reducing
// ordering will take it into account, otherwise it will be ignored.
bool ReorderProgramForSparseNormalCholesky(
SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
const ParameterBlockOrdering& parameter_block_ordering,
Program* program,
string* error);
} // namespace internal
} // namespace ceres
#endif // CERES_INTERNAL_REORDER_PROGRAM_

154
extern/libmv/third_party/ceres/internal/ceres/schur_complement_solver.cc vendored
View File

@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
@ -28,12 +28,13 @@
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#include "ceres/internal/port.h"
#include <algorithm>
#include <ctime>
#include <set>
#include <vector>
#include "Eigen/Dense"
#include "ceres/block_random_access_dense_matrix.h"
#include "ceres/block_random_access_matrix.h"
#include "ceres/block_random_access_sparse_matrix.h"
@ -42,7 +43,6 @@
#include "ceres/cxsparse.h"
#include "ceres/detect_structure.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/port.h"
#include "ceres/internal/scoped_ptr.h"
#include "ceres/lapack.h"
#include "ceres/linear_solver.h"
@ -51,6 +51,8 @@
#include "ceres/triplet_sparse_matrix.h"
#include "ceres/types.h"
#include "ceres/wall_time.h"
#include "Eigen/Dense"
#include "Eigen/SparseCore"
namespace ceres {
namespace internal {
@ -138,7 +140,8 @@ DenseSchurComplementSolver::SolveReducedLinearSystem(double* solution) {
.llt();
if (llt.info() != Eigen::Success) {
summary.termination_type = LINEAR_SOLVER_FAILURE;
summary.message = "Eigen LLT decomposition failed.";
summary.message =
"Eigen failure. Unable to perform dense Cholesky factorization.";
return summary;
}
@ -155,8 +158,6 @@ DenseSchurComplementSolver::SolveReducedLinearSystem(double* solution) {
return summary;
}
#if !defined(CERES_NO_SUITESPARSE) || !defined(CERES_NO_CXSPARE)
SparseSchurComplementSolver::SparseSchurComplementSolver(
const LinearSolver::Options& options)
: SchurComplementSolver(options),
@ -165,19 +166,15 @@ SparseSchurComplementSolver::SparseSchurComplementSolver(
}
SparseSchurComplementSolver::~SparseSchurComplementSolver() {
#ifndef CERES_NO_SUITESPARSE
if (factor_ != NULL) {
ss_.Free(factor_);
factor_ = NULL;
}
#endif // CERES_NO_SUITESPARSE
#ifndef CERES_NO_CXSPARSE
if (cxsparse_factor_ != NULL) {
cxsparse_.Free(cxsparse_factor_);
cxsparse_factor_ = NULL;
}
#endif // CERES_NO_CXSPARSE
}
// Determine the non-zero blocks in the Schur Complement matrix, and
@ -258,6 +255,8 @@ SparseSchurComplementSolver::SolveReducedLinearSystem(double* solution) {
return SolveReducedLinearSystemUsingSuiteSparse(solution);
case CX_SPARSE:
return SolveReducedLinearSystemUsingCXSparse(solution);
case EIGEN_SPARSE:
return SolveReducedLinearSystemUsingEigen(solution);
default:
LOG(FATAL) << "Unknown sparse linear algebra library : "
<< options().sparse_linear_algebra_library_type;
@ -266,13 +265,23 @@ SparseSchurComplementSolver::SolveReducedLinearSystem(double* solution) {
return LinearSolver::Summary();
}
#ifndef CERES_NO_SUITESPARSE
// Solve the system Sx = r, assuming that the matrix S is stored in a
// BlockRandomAccessSparseMatrix. The linear system is solved using
// CHOLMOD's sparse cholesky factorization routines.
LinearSolver::Summary
SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse(
double* solution) {
#ifdef CERES_NO_SUITESPARSE
LinearSolver::Summary summary;
summary.num_iterations = 0;
summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
summary.message = "Ceres was not built with SuiteSparse support. "
"Therefore, SPARSE_SCHUR cannot be used with SUITE_SPARSE";
return summary;
#else
LinearSolver::Summary summary;
summary.num_iterations = 0;
summary.termination_type = LINEAR_SOLVER_SUCCESS;
@ -326,6 +335,8 @@ SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse(
if (factor_ == NULL) {
ss_.Free(cholmod_lhs);
summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
// No need to set message as it has already been set by the
// symbolic analysis routines above.
return summary;
}
@ -335,6 +346,8 @@ SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse(
ss_.Free(cholmod_lhs);
if (summary.termination_type != LINEAR_SOLVER_SUCCESS) {
// No need to set message as it has already been set by the
// numeric factorization routine above.
return summary;
}
@ -346,6 +359,8 @@ SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse(
ss_.Free(cholmod_rhs);
if (cholmod_solution == NULL) {
summary.message =
"SuiteSparse failure. Unable to perform triangular solve.";
summary.termination_type = LINEAR_SOLVER_FAILURE;
return summary;
}
@ -354,23 +369,26 @@ SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse(
= VectorRef(static_cast<double*>(cholmod_solution->x), num_rows);
ss_.Free(cholmod_solution);
return summary;
}
#else
LinearSolver::Summary
SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse(
double* solution) {
LOG(FATAL) << "No SuiteSparse support in Ceres.";
return LinearSolver::Summary();
}
#endif // CERES_NO_SUITESPARSE
}
#ifndef CERES_NO_CXSPARSE
// Solve the system Sx = r, assuming that the matrix S is stored in a
// BlockRandomAccessSparseMatrix. The linear system is solved using
// CXSparse's sparse cholesky factorization routines.
LinearSolver::Summary
SparseSchurComplementSolver::SolveReducedLinearSystemUsingCXSparse(
double* solution) {
#ifdef CERES_NO_CXSPARSE
LinearSolver::Summary summary;
summary.num_iterations = 0;
summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
summary.message = "Ceres was not built with CXSparse support. "
"Therefore, SPARSE_SCHUR cannot be used with CX_SPARSE";
return summary;
#else
LinearSolver::Summary summary;
summary.num_iterations = 0;
summary.termination_type = LINEAR_SOLVER_SUCCESS;
@ -407,16 +425,94 @@ SparseSchurComplementSolver::SolveReducedLinearSystemUsingCXSparse(
cxsparse_.Free(lhs);
return summary;
}
#else
LinearSolver::Summary
SparseSchurComplementSolver::SolveReducedLinearSystemUsingCXSparse(
double* solution) {
LOG(FATAL) << "No CXSparse support in Ceres.";
return LinearSolver::Summary();
}
#endif // CERES_NO_CXPARSE
}
// Solve the system Sx = r, assuming that the matrix S is stored in a
// BlockRandomAccessSparseMatrix. The linear system is solved using
// Eigen's sparse cholesky factorization routines.
LinearSolver::Summary
SparseSchurComplementSolver::SolveReducedLinearSystemUsingEigen(
double* solution) {
#ifndef CERES_USE_EIGEN_SPARSE
LinearSolver::Summary summary;
summary.num_iterations = 0;
summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
summary.message =
"SPARSE_SCHUR cannot be used with EIGEN_SPARSE. "
"Ceres was not built with support for "
"Eigen's SimplicialLDLT decomposition. "
"This requires enabling building with -DEIGENSPARSE=ON.";
return summary;
#else
EventLogger event_logger("SchurComplementSolver::EigenSolve");
LinearSolver::Summary summary;
summary.num_iterations = 0;
summary.termination_type = LINEAR_SOLVER_SUCCESS;
summary.message = "Success.";
// Extract the TripletSparseMatrix that is used for actually storing S.
TripletSparseMatrix* tsm =
const_cast<TripletSparseMatrix*>(
down_cast<const BlockRandomAccessSparseMatrix*>(lhs())->matrix());
const int num_rows = tsm->num_rows();
// The case where there are no f blocks, and the system is block
// diagonal.
if (num_rows == 0) {
return summary;
}
// This is an upper triangular matrix.
CompressedRowSparseMatrix crsm(*tsm);
// Map this to a column major, lower triangular matrix.
Eigen::MappedSparseMatrix<double, Eigen::ColMajor> eigen_lhs(
crsm.num_rows(),
crsm.num_rows(),
crsm.num_nonzeros(),
crsm.mutable_rows(),
crsm.mutable_cols(),
crsm.mutable_values());
event_logger.AddEvent("ToCompressedRowSparseMatrix");
// Compute symbolic factorization if one does not exist.
if (simplicial_ldlt_.get() == NULL) {
simplicial_ldlt_.reset(new SimplicialLDLT);
// This ordering is quite bad. The scalar ordering produced by the
// AMD algorithm is quite bad and can be an order of magnitude
// worse than the one computed using the block version of the
// algorithm.
simplicial_ldlt_->analyzePattern(eigen_lhs);
event_logger.AddEvent("Analysis");
if (simplicial_ldlt_->info() != Eigen::Success) {
summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
summary.message =
"Eigen failure. Unable to find symbolic factorization.";
return summary;
}
}
simplicial_ldlt_->factorize(eigen_lhs);
event_logger.AddEvent("Factorize");
if (simplicial_ldlt_->info() != Eigen::Success) {
summary.termination_type = LINEAR_SOLVER_FAILURE;
summary.message = "Eigen failure. Unable to find numeric factoriztion.";
return summary;
}
VectorRef(solution, num_rows) =
simplicial_ldlt_->solve(ConstVectorRef(rhs(), num_rows));
event_logger.AddEvent("Solve");
if (simplicial_ldlt_->info() != Eigen::Success) {
summary.termination_type = LINEAR_SOLVER_FAILURE;
summary.message = "Eigen failure. Unable to do triangular solve.";
}
return summary;
#endif // CERES_USE_EIGEN_SPARSE
}
#endif // !defined(CERES_NO_SUITESPARSE) || !defined(CERES_NO_CXSPARE)
} // namespace internal
} // namespace ceres

29
extern/libmv/third_party/ceres/internal/ceres/schur_complement_solver.h vendored
View File

@ -35,6 +35,8 @@
#include <utility>
#include <vector>
#include "ceres/internal/port.h"
#include "ceres/block_random_access_matrix.h"
#include "ceres/block_sparse_matrix.h"
#include "ceres/block_structure.h"
@ -45,6 +47,11 @@
#include "ceres/internal/scoped_ptr.h"
#include "ceres/types.h"
#ifdef CERES_USE_EIGEN_SPARSE
#include "Eigen/SparseCholesky"
#include "Eigen/OrderingMethods"
#endif
namespace ceres {
namespace internal {
@ -153,7 +160,6 @@ class DenseSchurComplementSolver : public SchurComplementSolver {
CERES_DISALLOW_COPY_AND_ASSIGN(DenseSchurComplementSolver);
};
#if !defined(CERES_NO_SUITESPARSE) || !defined(CERES_NO_CXSPARE)
// Sparse Cholesky factorization based solver.
class SparseSchurComplementSolver : public SchurComplementSolver {
public:
@ -168,6 +174,8 @@ class SparseSchurComplementSolver : public SchurComplementSolver {
double* solution);
LinearSolver::Summary SolveReducedLinearSystemUsingCXSparse(
double* solution);
LinearSolver::Summary SolveReducedLinearSystemUsingEigen(
double* solution);
// Size of the blocks in the Schur complement.
vector<int> blocks_;
@ -180,10 +188,27 @@ class SparseSchurComplementSolver : public SchurComplementSolver {
CXSparse cxsparse_;
// Cached factorization
cs_dis* cxsparse_factor_;
#ifdef CERES_USE_EIGEN_SPARSE
// The preprocessor gymnastics here are dealing with the fact that
// before version 3.2.2, Eigen did not support a third template
// parameter to specify the ordering.
#if EIGEN_VERSION_AT_LEAST(3,2,2)
typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>, Eigen::Lower,
Eigen::NaturalOrdering<int> >
SimplicialLDLT;
#else
typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>, Eigen::Lower>
SimplicialLDLT;
#endif
scoped_ptr<SimplicialLDLT> simplicial_ldlt_;
#endif
CERES_DISALLOW_COPY_AND_ASSIGN(SparseSchurComplementSolver);
};
#endif // !defined(CERES_NO_SUITESPARSE) || !defined(CERES_NO_CXSPARE)
} // namespace internal
} // namespace ceres

2
extern/libmv/third_party/ceres/internal/ceres/single_linkage_clustering.cc vendored
View File

@ -44,7 +44,7 @@ namespace internal {
int ComputeSingleLinkageClustering(
const SingleLinkageClusteringOptions& options,
const Graph<int>& graph,
const WeightedGraph<int>& graph,
HashMap<int, int>* membership) {
CHECK_NOTNULL(membership)->clear();

2
extern/libmv/third_party/ceres/internal/ceres/single_linkage_clustering.h vendored
View File

@ -64,7 +64,7 @@ struct SingleLinkageClusteringOptions {
// identified by the algorithm.
int ComputeSingleLinkageClustering(
const SingleLinkageClusteringOptions& options,
const Graph<int>& graph,
const WeightedGraph<int>& graph,
HashMap<int, int>* membership);
} // namespace internal

38
extern/libmv/third_party/ceres/internal/ceres/small_blas.h vendored
View File

@ -42,30 +42,6 @@
namespace ceres {
namespace internal {
// Remove the ".noalias()" annotation from the matrix matrix
// mutliplies to produce a correct build with the Android NDK,
// including versions 6, 7, 8, and 8b, when built with STLPort and the
// non-standalone toolchain (i.e. ndk-build). This appears to be a
// compiler bug; if the workaround is not in place, the line
//
// block.noalias() -= A * B;
//
// gets compiled to
//
// block.noalias() += A * B;
//
// which breaks schur elimination. Introducing a temporary by removing the
// .noalias() annotation causes the issue to disappear. Tracking this
// issue down was tricky, since the test suite doesn't run when built with
// the non-standalone toolchain.
//
// TODO(keir): Make a reproduction case for this and send it upstream.
#ifdef CERES_WORK_AROUND_ANDROID_NDK_COMPILER_BUG
#define CERES_MAYBE_NOALIAS
#else
#define CERES_MAYBE_NOALIAS .noalias()
#endif
// The following three macros are used to share code and reduce
// template junk across the various GEMM variants.
#define CERES_GEMM_BEGIN(name) \
@ -168,11 +144,11 @@ CERES_GEMM_BEGIN(MatrixMatrixMultiplyEigen) {
block(Cref, start_row_c, start_col_c, num_row_a, num_col_b);
if (kOperation > 0) {
block CERES_MAYBE_NOALIAS += Aref * Bref;
block.noalias() += Aref * Bref;
} else if (kOperation < 0) {
block CERES_MAYBE_NOALIAS -= Aref * Bref;
block.noalias() -= Aref * Bref;
} else {
block CERES_MAYBE_NOALIAS = Aref * Bref;
block.noalias() = Aref * Bref;
}
}
@ -228,11 +204,11 @@ CERES_GEMM_BEGIN(MatrixTransposeMatrixMultiplyEigen) {
start_row_c, start_col_c,
num_col_a, num_col_b);
if (kOperation > 0) {
block CERES_MAYBE_NOALIAS += Aref.transpose() * Bref;
block.noalias() += Aref.transpose() * Bref;
} else if (kOperation < 0) {
block CERES_MAYBE_NOALIAS -= Aref.transpose() * Bref;
block.noalias() -= Aref.transpose() * Bref;
} else {
block CERES_MAYBE_NOALIAS = Aref.transpose() * Bref;
block.noalias() = Aref.transpose() * Bref;
}
}
@ -394,8 +370,6 @@ inline void MatrixTransposeVectorMultiply(const double* A,
#endif // CERES_NO_CUSTOM_BLAS
}
#undef CERES_MAYBE_NOALIAS
#undef CERES_GEMM_BEGIN
#undef CERES_GEMM_EIGEN_HEADER
#undef CERES_GEMM_NAIVE_HEADER

491
extern/libmv/third_party/ceres/internal/ceres/solver.cc vendored
View File

@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
@ -31,17 +31,274 @@
#include "ceres/solver.h"
#include <algorithm>
#include <sstream> // NOLINT
#include <vector>
#include "ceres/gradient_checking_cost_function.h"
#include "ceres/internal/port.h"
#include "ceres/parameter_block_ordering.h"
#include "ceres/preprocessor.h"
#include "ceres/problem.h"
#include "ceres/problem_impl.h"
#include "ceres/program.h"
#include "ceres/solver_impl.h"
#include "ceres/solver_utils.h"
#include "ceres/stringprintf.h"
#include "ceres/types.h"
#include "ceres/wall_time.h"
namespace ceres {
namespace {
#define OPTION_OP(x, y, OP) \
if (!(options.x OP y)) { \
std::stringstream ss; \
ss << "Invalid configuration. "; \
ss << string("Solver::Options::" #x " = ") << options.x << ". "; \
ss << "Violated constraint: "; \
ss << string("Solver::Options::" #x " " #OP " "#y); \
*error = ss.str(); \
return false; \
}
#define OPTION_OP_OPTION(x, y, OP) \
if (!(options.x OP options.y)) { \
std::stringstream ss; \
ss << "Invalid configuration. "; \
ss << string("Solver::Options::" #x " = ") << options.x << ". "; \
ss << string("Solver::Options::" #y " = ") << options.y << ". "; \
ss << "Violated constraint: "; \
ss << string("Solver::Options::" #x); \
ss << string(#OP " Solver::Options::" #y "."); \
*error = ss.str(); \
return false; \
}
#define OPTION_GE(x, y) OPTION_OP(x, y, >=);
#define OPTION_GT(x, y) OPTION_OP(x, y, >);
#define OPTION_LE(x, y) OPTION_OP(x, y, <=);
#define OPTION_LT(x, y) OPTION_OP(x, y, <);
#define OPTION_LE_OPTION(x, y) OPTION_OP_OPTION(x, y, <=)
#define OPTION_LT_OPTION(x, y) OPTION_OP_OPTION(x, y, <)
bool CommonOptionsAreValid(const Solver::Options& options, string* error) {
OPTION_GE(max_num_iterations, 0);
OPTION_GE(max_solver_time_in_seconds, 0.0);
OPTION_GE(function_tolerance, 0.0);
OPTION_GE(gradient_tolerance, 0.0);
OPTION_GE(parameter_tolerance, 0.0);
OPTION_GT(num_threads, 0);
OPTION_GT(num_linear_solver_threads, 0);
if (options.check_gradients) {
OPTION_GT(gradient_check_relative_precision, 0.0);
OPTION_GT(numeric_derivative_relative_step_size, 0.0);
}
return true;
}
bool TrustRegionOptionsAreValid(const Solver::Options& options, string* error) {
OPTION_GT(initial_trust_region_radius, 0.0);
OPTION_GT(min_trust_region_radius, 0.0);
OPTION_GT(max_trust_region_radius, 0.0);
OPTION_LE_OPTION(min_trust_region_radius, max_trust_region_radius);
OPTION_LE_OPTION(min_trust_region_radius, initial_trust_region_radius);
OPTION_LE_OPTION(initial_trust_region_radius, max_trust_region_radius);
OPTION_GE(min_relative_decrease, 0.0);
OPTION_GE(min_lm_diagonal, 0.0);
OPTION_GE(max_lm_diagonal, 0.0);
OPTION_LE_OPTION(min_lm_diagonal, max_lm_diagonal);
OPTION_GE(max_num_consecutive_invalid_steps, 0);
OPTION_GT(eta, 0.0);
OPTION_GE(min_linear_solver_iterations, 0);
OPTION_GE(max_linear_solver_iterations, 1);
OPTION_LE_OPTION(min_linear_solver_iterations, max_linear_solver_iterations);
if (options.use_inner_iterations) {
OPTION_GE(inner_iteration_tolerance, 0.0);
}
if (options.use_nonmonotonic_steps) {
OPTION_GT(max_consecutive_nonmonotonic_steps, 0);
}
if (options.preconditioner_type == CLUSTER_JACOBI &&
options.sparse_linear_algebra_library_type != SUITE_SPARSE) {
*error = "CLUSTER_JACOBI requires "
"Solver::Options::sparse_linear_algebra_library_type to be "
"SUITE_SPARSE";
return false;
}
if (options.preconditioner_type == CLUSTER_TRIDIAGONAL &&
options.sparse_linear_algebra_library_type != SUITE_SPARSE) {
*error = "CLUSTER_TRIDIAGONAL requires "
"Solver::Options::sparse_linear_algebra_library_type to be "
"SUITE_SPARSE";
return false;
}
#ifdef CERES_NO_LAPACK
if (options.dense_linear_algebra_library_type == LAPACK) {
if (options.linear_solver_type == DENSE_NORMAL_CHOLESKY) {
*error = "Can't use DENSE_NORMAL_CHOLESKY with LAPACK because "
"LAPACK was not enabled when Ceres was built.";
return false;
}
if (options.linear_solver_type == DENSE_QR) {
*error = "Can't use DENSE_QR with LAPACK because "
"LAPACK was not enabled when Ceres was built.";
return false;
}
if (options.linear_solver_type == DENSE_SCHUR) {
*error = "Can't use DENSE_SCHUR with LAPACK because "
"LAPACK was not enabled when Ceres was built.";
return false;
}
}
#endif
#ifdef CERES_NO_SUITESPARSE
if (options.sparse_linear_algebra_library_type == SUITE_SPARSE) {
if (options.linear_solver_type == SPARSE_NORMAL_CHOLESKY) {
*error = "Can't use SPARSE_NORMAL_CHOLESKY with SUITESPARSE because "
"SuiteSparse was not enabled when Ceres was built.";
return false;
}
if (options.linear_solver_type == SPARSE_SCHUR) {
*error = "Can't use SPARSE_SCHUR with SUITESPARSE because "
"SuiteSparse was not enabled when Ceres was built.";
return false;
}
if (options.preconditioner_type == CLUSTER_JACOBI) {
*error = "CLUSTER_JACOBI preconditioner not supported. "
"SuiteSparse was not enabled when Ceres was built.";
return false;
}
if (options.preconditioner_type == CLUSTER_TRIDIAGONAL) {
*error = "CLUSTER_TRIDIAGONAL preconditioner not supported. "
"SuiteSparse was not enabled when Ceres was built.";
return false;
}
}
#endif
#ifdef CERES_NO_CXSPARSE
if (options.sparse_linear_algebra_library_type == CX_SPARSE) {
if (options.linear_solver_type == SPARSE_NORMAL_CHOLESKY) {
*error = "Can't use SPARSE_NORMAL_CHOLESKY with CX_SPARSE because "
"CXSparse was not enabled when Ceres was built.";
return false;
}
if (options.linear_solver_type == SPARSE_SCHUR) {
*error = "Can't use SPARSE_SCHUR with CX_SPARSE because "
"CXSparse was not enabled when Ceres was built.";
return false;
}
}
#endif
#ifndef CERES_USE_EIGEN_SPARSE
if (options.sparse_linear_algebra_library_type == EIGEN_SPARSE) {
if (options.linear_solver_type == SPARSE_NORMAL_CHOLESKY) {
*error = "Can't use SPARSE_NORMAL_CHOLESKY with EIGEN_SPARSE because "
"Eigen's sparse linear algebra was not enabled when Ceres was "
"built.";
return false;
}
if (options.linear_solver_type == SPARSE_SCHUR) {
*error = "Can't use SPARSE_SCHUR with EIGEN_SPARSE because "
"Eigen's sparse linear algebra was not enabled when Ceres was "
"built.";
return false;
}
}
#endif
if (options.trust_region_strategy_type == DOGLEG) {
if (options.linear_solver_type == ITERATIVE_SCHUR ||
options.linear_solver_type == CGNR) {
*error = "DOGLEG only supports exact factorization based linear "
"solvers. If you want to use an iterative solver please "
"use LEVENBERG_MARQUARDT as the trust_region_strategy_type";
return false;
}
}
if (options.trust_region_minimizer_iterations_to_dump.size() > 0 &&
options.trust_region_problem_dump_format_type != CONSOLE &&
options.trust_region_problem_dump_directory.empty()) {
*error = "Solver::Options::trust_region_problem_dump_directory is empty.";
return false;
}
if (options.dynamic_sparsity &&
options.linear_solver_type != SPARSE_NORMAL_CHOLESKY) {
*error = "Dynamic sparsity is only supported with SPARSE_NORMAL_CHOLESKY.";
return false;
}
return true;
}
bool LineSearchOptionsAreValid(const Solver::Options& options, string* error) {
OPTION_GT(max_lbfgs_rank, 0);
OPTION_GT(min_line_search_step_size, 0.0);
OPTION_GT(max_line_search_step_contraction, 0.0);
OPTION_LT(max_line_search_step_contraction, 1.0);
OPTION_LT_OPTION(max_line_search_step_contraction,
min_line_search_step_contraction);
OPTION_LE(min_line_search_step_contraction, 1.0);
OPTION_GT(max_num_line_search_step_size_iterations, 0);
OPTION_GT(line_search_sufficient_function_decrease, 0.0);
OPTION_LT_OPTION(line_search_sufficient_function_decrease,
line_search_sufficient_curvature_decrease);
OPTION_LT(line_search_sufficient_curvature_decrease, 1.0);
OPTION_GT(max_line_search_step_expansion, 1.0);
if ((options.line_search_direction_type == ceres::BFGS ||
options.line_search_direction_type == ceres::LBFGS) &&
options.line_search_type != ceres::WOLFE) {
*error =
string("Invalid configuration: Solver::Options::line_search_type = ")
+ string(LineSearchTypeToString(options.line_search_type))
+ string(". When using (L)BFGS, "
"Solver::Options::line_search_type must be set to WOLFE.");
return false;
}
// Warn user if they have requested BISECTION interpolation, but constraints
// on max/min step size change during line search prevent bisection scaling
// from occurring. Warn only, as this is likely a user mistake, but one which
// does not prevent us from continuing.
LOG_IF(WARNING,
(options.line_search_interpolation_type == ceres::BISECTION &&
(options.max_line_search_step_contraction > 0.5 ||
options.min_line_search_step_contraction < 0.5)))
<< "Line search interpolation type is BISECTION, but specified "
<< "max_line_search_step_contraction: "
<< options.max_line_search_step_contraction << ", and "
<< "min_line_search_step_contraction: "
<< options.min_line_search_step_contraction
<< ", prevent bisection (0.5) scaling, continuing with solve regardless.";
return true;
}
#undef OPTION_OP
#undef OPTION_OP_OPTION
#undef OPTION_GT
#undef OPTION_GE
#undef OPTION_LE
#undef OPTION_LT
#undef OPTION_LE_OPTION
#undef OPTION_LT_OPTION
void StringifyOrdering(const vector<int>& ordering, string* report) {
if (ordering.size() == 0) {
internal::StringAppendF(report, "AUTOMATIC");
@ -54,19 +311,211 @@ void StringifyOrdering(const vector<int>& ordering, string* report) {
internal::StringAppendF(report, "%d", ordering.back());
}
void SummarizeGivenProgram(const internal::Program& program,
Solver::Summary* summary) {
summary->num_parameter_blocks = program.NumParameterBlocks();
summary->num_parameters = program.NumParameters();
summary->num_effective_parameters = program.NumEffectiveParameters();
summary->num_residual_blocks = program.NumResidualBlocks();
summary->num_residuals = program.NumResiduals();
}
void SummarizeReducedProgram(const internal::Program& program,
Solver::Summary* summary) {
summary->num_parameter_blocks_reduced = program.NumParameterBlocks();
summary->num_parameters_reduced = program.NumParameters();
summary->num_effective_parameters_reduced = program.NumEffectiveParameters();
summary->num_residual_blocks_reduced = program.NumResidualBlocks();
summary->num_residuals_reduced = program.NumResiduals();
}
void PreSolveSummarize(const Solver::Options& options,
const internal::ProblemImpl* problem,
Solver::Summary* summary) {
SummarizeGivenProgram(problem->program(), summary);
internal::OrderingToGroupSizes(options.linear_solver_ordering.get(),
&(summary->linear_solver_ordering_given));
internal::OrderingToGroupSizes(options.inner_iteration_ordering.get(),
&(summary->inner_iteration_ordering_given));
summary->dense_linear_algebra_library_type = options.dense_linear_algebra_library_type; // NOLINT
summary->dogleg_type = options.dogleg_type;
summary->inner_iteration_time_in_seconds = 0.0;
summary->inner_iterations_given = options.use_inner_iterations;
summary->line_search_direction_type = options.line_search_direction_type; // NOLINT
summary->line_search_interpolation_type = options.line_search_interpolation_type; // NOLINT
summary->line_search_type = options.line_search_type;
summary->linear_solver_type_given = options.linear_solver_type;
summary->max_lbfgs_rank = options.max_lbfgs_rank;
summary->minimizer_type = options.minimizer_type;
summary->nonlinear_conjugate_gradient_type = options.nonlinear_conjugate_gradient_type; // NOLINT
summary->num_linear_solver_threads_given = options.num_linear_solver_threads; // NOLINT
summary->num_threads_given = options.num_threads;
summary->preconditioner_type_given = options.preconditioner_type;
summary->sparse_linear_algebra_library_type = options.sparse_linear_algebra_library_type; // NOLINT
summary->trust_region_strategy_type = options.trust_region_strategy_type; // NOLINT
summary->visibility_clustering_type = options.visibility_clustering_type; // NOLINT
}
void PostSolveSummarize(const internal::PreprocessedProblem& pp,
Solver::Summary* summary) {
internal::OrderingToGroupSizes(pp.options.linear_solver_ordering.get(),
&(summary->linear_solver_ordering_used));
internal::OrderingToGroupSizes(pp.options.inner_iteration_ordering.get(),
&(summary->inner_iteration_ordering_used));
summary->inner_iterations_used = pp.inner_iteration_minimizer.get() != NULL; // NOLINT
summary->linear_solver_type_used = pp.options.linear_solver_type;
summary->num_linear_solver_threads_used = pp.options.num_linear_solver_threads; // NOLINT
summary->num_threads_used = pp.options.num_threads;
summary->preconditioner_type_used = pp.options.preconditioner_type; // NOLINT
internal::SetSummaryFinalCost(summary);
if (pp.reduced_program.get() != NULL) {
SummarizeReducedProgram(*pp.reduced_program, summary);
}
// It is possible that no evaluator was created. This would be the
// case if the preprocessor failed, or if the reduced problem did
// not contain any parameter blocks. Thus, only extract the
// evaluator statistics if one exists.
if (pp.evaluator.get() != NULL) {
const map<string, double>& evaluator_time_statistics =
pp.evaluator->TimeStatistics();
summary->residual_evaluation_time_in_seconds =
FindWithDefault(evaluator_time_statistics, "Evaluator::Residual", 0.0);
summary->jacobian_evaluation_time_in_seconds =
FindWithDefault(evaluator_time_statistics, "Evaluator::Jacobian", 0.0);
}
// Again, like the evaluator, there may or may not be a linear
// solver from which we can extract run time statistics. In
// particular the line search solver does not use a linear solver.
if (pp.linear_solver.get() != NULL) {
const map<string, double>& linear_solver_time_statistics =
pp.linear_solver->TimeStatistics();
summary->linear_solver_time_in_seconds =
FindWithDefault(linear_solver_time_statistics,
"LinearSolver::Solve",
0.0);
}
}
void Minimize(internal::PreprocessedProblem* pp,
Solver::Summary* summary) {
using internal::Program;
using internal::scoped_ptr;
using internal::Minimizer;
Program* program = pp->reduced_program.get();
if (pp->reduced_program->NumParameterBlocks() == 0) {
summary->message = "Function tolerance reached. "
"No non-constant parameter blocks found.";
summary->termination_type = CONVERGENCE;
VLOG_IF(1, pp->options.logging_type != SILENT) << summary->message;
summary->initial_cost = summary->fixed_cost;
summary->final_cost = summary->fixed_cost;
return;
}
scoped_ptr<Minimizer> minimizer(
Minimizer::Create(pp->options.minimizer_type));
minimizer->Minimize(pp->minimizer_options,
pp->reduced_parameters.data(),
summary);
if (summary->IsSolutionUsable()) {
program->StateVectorToParameterBlocks(pp->reduced_parameters.data());
program->CopyParameterBlockStateToUserState();
}
}
} // namespace
bool Solver::Options::IsValid(string* error) const {
if (!CommonOptionsAreValid(*this, error)) {
return false;
}
if (minimizer_type == TRUST_REGION) {
return TrustRegionOptionsAreValid(*this, error);
}
CHECK_EQ(minimizer_type, LINE_SEARCH);
return LineSearchOptionsAreValid(*this, error);
}
Solver::~Solver() {}
void Solver::Solve(const Solver::Options& options,
Problem* problem,
Solver::Summary* summary) {
double start_time_seconds = internal::WallTimeInSeconds();
internal::ProblemImpl* problem_impl =
CHECK_NOTNULL(problem)->problem_impl_.get();
internal::SolverImpl::Solve(options, problem_impl, summary);
summary->total_time_in_seconds =
internal::WallTimeInSeconds() - start_time_seconds;
using internal::PreprocessedProblem;
using internal::Preprocessor;
using internal::ProblemImpl;
using internal::Program;
using internal::scoped_ptr;
using internal::WallTimeInSeconds;
CHECK_NOTNULL(problem);
CHECK_NOTNULL(summary);
double start_time = WallTimeInSeconds();
*summary = Summary();
if (!options.IsValid(&summary->message)) {
LOG(ERROR) << "Terminating: " << summary->message;
return;
}
ProblemImpl* problem_impl = problem->problem_impl_.get();
Program* program = problem_impl->mutable_program();
PreSolveSummarize(options, problem_impl, summary);
// Make sure that all the parameter blocks states are set to the
// values provided by the user.
program->SetParameterBlockStatePtrsToUserStatePtrs();
scoped_ptr<internal::ProblemImpl> gradient_checking_problem;
if (options.check_gradients) {
gradient_checking_problem.reset(
CreateGradientCheckingProblemImpl(
problem_impl,
options.numeric_derivative_relative_step_size,
options.gradient_check_relative_precision));
problem_impl = gradient_checking_problem.get();
program = problem_impl->mutable_program();
}
scoped_ptr<Preprocessor> preprocessor(
Preprocessor::Create(options.minimizer_type));
PreprocessedProblem pp;
const bool status = preprocessor->Preprocess(options, problem_impl, &pp);
summary->fixed_cost = pp.fixed_cost;
summary->preprocessor_time_in_seconds = WallTimeInSeconds() - start_time;
if (status) {
const double minimizer_start_time = WallTimeInSeconds();
Minimize(&pp, summary);
summary->minimizer_time_in_seconds =
WallTimeInSeconds() - minimizer_start_time;
} else {
summary->message = pp.error;
}
const double postprocessor_start_time = WallTimeInSeconds();
problem_impl = problem->problem_impl_.get();
program = problem_impl->mutable_program();
// On exit, ensure that the parameter blocks again point at the user
// provided values and the parameter blocks are numbered according
// to their position in the original user provided program.
program->SetParameterBlockStatePtrsToUserStatePtrs();
program->SetParameterOffsetsAndIndex();
PostSolveSummarize(pp, summary);
summary->postprocessor_time_in_seconds =
WallTimeInSeconds() - postprocessor_start_time;
summary->total_time_in_seconds = WallTimeInSeconds() - start_time;
}
void Solve(const Solver::Options& options,
@ -114,7 +563,8 @@ Solver::Summary::Summary()
linear_solver_type_used(SPARSE_NORMAL_CHOLESKY),
inner_iterations_given(false),
inner_iterations_used(false),
preconditioner_type(IDENTITY),
preconditioner_type_given(IDENTITY),
preconditioner_type_used(IDENTITY),
visibility_clustering_type(CANONICAL_VIEWS),
trust_region_strategy_type(LEVENBERG_MARQUARDT),
dense_linear_algebra_library_type(EIGEN),
@ -142,10 +592,9 @@ string Solver::Summary::BriefReport() const {
};
string Solver::Summary::FullReport() const {
string report =
"\n"
"Ceres Solver Report\n"
"-------------------\n";
using internal::VersionString;
string report = string("\nSolver Summary (v " + VersionString() + ")\n\n");
StringAppendF(&report, "%45s %21s\n", "Original", "Reduced");
StringAppendF(&report, "Parameter blocks % 25d% 25d\n",
@ -177,8 +626,8 @@ string Solver::Summary::FullReport() const {
if (linear_solver_type_used == SPARSE_NORMAL_CHOLESKY ||
linear_solver_type_used == SPARSE_SCHUR ||
(linear_solver_type_used == ITERATIVE_SCHUR &&
(preconditioner_type == CLUSTER_JACOBI ||
preconditioner_type == CLUSTER_TRIDIAGONAL))) {
(preconditioner_type_used == CLUSTER_JACOBI ||
preconditioner_type_used == CLUSTER_TRIDIAGONAL))) {
StringAppendF(&report, "\nSparse linear algebra library %15s\n",
SparseLinearAlgebraLibraryTypeToString(
sparse_linear_algebra_library_type));
@ -205,12 +654,12 @@ string Solver::Summary::FullReport() const {
if (linear_solver_type_given == CGNR ||
linear_solver_type_given == ITERATIVE_SCHUR) {
StringAppendF(&report, "Preconditioner %25s%25s\n",
PreconditionerTypeToString(preconditioner_type),
PreconditionerTypeToString(preconditioner_type));
PreconditionerTypeToString(preconditioner_type_given),
PreconditionerTypeToString(preconditioner_type_used));
}
if (preconditioner_type == CLUSTER_JACOBI ||
preconditioner_type == CLUSTER_TRIDIAGONAL) {
if (preconditioner_type_used == CLUSTER_JACOBI ||
preconditioner_type_used == CLUSTER_TRIDIAGONAL) {
StringAppendF(&report, "Visibility clustering%24s%25s\n",
VisibilityClusteringTypeToString(
visibility_clustering_type),
@ -345,9 +794,7 @@ string Solver::Summary::FullReport() const {
};
bool Solver::Summary::IsSolutionUsable() const {
return (termination_type == CONVERGENCE ||
termination_type == NO_CONVERGENCE ||
termination_type == USER_SUCCESS);
return internal::IsSolutionUsable(*this);
}
} // namespace ceres

1851
extern/libmv/third_party/ceres/internal/ceres/solver_impl.cc vendored
View File

File diff suppressed because it is too large Load Diff

228
extern/libmv/third_party/ceres/internal/ceres/solver_impl.h vendored
View File

@ -1,228 +0,0 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: keir@google.com (Keir Mierle)
#ifndef CERES_INTERNAL_SOLVER_IMPL_H_
#define CERES_INTERNAL_SOLVER_IMPL_H_
#include <set>
#include <string>
#include <vector>
#include "ceres/internal/port.h"
#include "ceres/ordered_groups.h"
#include "ceres/problem_impl.h"
#include "ceres/solver.h"
namespace ceres {
namespace internal {
class CoordinateDescentMinimizer;
class Evaluator;
class LinearSolver;
class Program;
class TripletSparseMatrix;
class SolverImpl {
public:
// Mirrors the interface in solver.h, but exposes implementation
// details for testing internally.
static void Solve(const Solver::Options& options,
ProblemImpl* problem_impl,
Solver::Summary* summary);
static void TrustRegionSolve(const Solver::Options& options,
ProblemImpl* problem_impl,
Solver::Summary* summary);
// Run the TrustRegionMinimizer for the given evaluator and configuration.
static void TrustRegionMinimize(
const Solver::Options &options,
Program* program,
CoordinateDescentMinimizer* inner_iteration_minimizer,
Evaluator* evaluator,
LinearSolver* linear_solver,
Solver::Summary* summary);
static void LineSearchSolve(const Solver::Options& options,
ProblemImpl* problem_impl,
Solver::Summary* summary);
// Run the LineSearchMinimizer for the given evaluator and configuration.
static void LineSearchMinimize(const Solver::Options &options,
Program* program,
Evaluator* evaluator,
Solver::Summary* summary);
// Create the transformed Program, which has all the fixed blocks
// and residuals eliminated, and in the case of automatic schur
// ordering, has the E blocks first in the resulting program, with
// options.num_eliminate_blocks set appropriately.
//
// If fixed_cost is not NULL, the residual blocks that are removed
// are evaluated and the sum of their cost is returned in fixed_cost.
static Program* CreateReducedProgram(Solver::Options* options,
ProblemImpl* problem_impl,
double* fixed_cost,
string* message);
// Create the appropriate linear solver, taking into account any
// config changes decided by CreateTransformedProgram(). The
// selected linear solver, which may be different from what the user
// selected; consider the case that the remaining elimininated
// blocks is zero after removing fixed blocks.
static LinearSolver* CreateLinearSolver(Solver::Options* options,
string* message);
// Reorder the residuals for program, if necessary, so that the
// residuals involving e block (i.e., the first num_eliminate_block
// parameter blocks) occur together. This is a necessary condition
// for the Schur eliminator.
static bool LexicographicallyOrderResidualBlocks(
const int num_eliminate_blocks,
Program* program,
string* message);
// Create the appropriate evaluator for the transformed program.
static Evaluator* CreateEvaluator(
const Solver::Options& options,
const ProblemImpl::ParameterMap& parameter_map,
Program* program,
string* message);
// Remove the fixed or unused parameter blocks and residuals
// depending only on fixed parameters from the program.
//
// If either linear_solver_ordering or inner_iteration_ordering are
// not NULL, the constant parameter blocks are removed from them
// too.
//
// If fixed_cost is not NULL, the residual blocks that are removed
// are evaluated and the sum of their cost is returned in
// fixed_cost.
//
// If a failure is encountered, the function returns false with a
// description of the failure in message.
static bool RemoveFixedBlocksFromProgram(
Program* program,
ParameterBlockOrdering* linear_solver_ordering,
ParameterBlockOrdering* inner_iteration_ordering,
double* fixed_cost,
string* message);
static bool IsOrderingValid(const Solver::Options& options,
const ProblemImpl* problem_impl,
string* message);
static bool IsParameterBlockSetIndependent(
const set<double*>& parameter_block_ptrs,
const vector<ResidualBlock*>& residual_blocks);
static CoordinateDescentMinimizer* CreateInnerIterationMinimizer(
const Solver::Options& options,
const Program& program,
const ProblemImpl::ParameterMap& parameter_map,
Solver::Summary* summary);
// If the linear solver is of Schur type, then replace it with the
// closest equivalent linear solver. This is done when the user
// requested a Schur type solver but the problem structure makes it
// impossible to use one.
//
// If the linear solver is not of Schur type, the function is a
// no-op.
static void AlternateLinearSolverForSchurTypeLinearSolver(
Solver::Options* options);
// Create a TripletSparseMatrix which contains the zero-one
// structure corresponding to the block sparsity of the transpose of
// the Jacobian matrix.
//
// Caller owns the result.
static TripletSparseMatrix* CreateJacobianBlockSparsityTranspose(
const Program* program);
// Reorder the parameter blocks in program using the ordering
static bool ApplyUserOrdering(
const ProblemImpl::ParameterMap& parameter_map,
const ParameterBlockOrdering* parameter_block_ordering,
Program* program,
string* message);
// Sparse cholesky factorization routines when doing the sparse
// cholesky factorization of the Jacobian matrix, reorders its
// columns to reduce the fill-in. Compute this permutation and
// re-order the parameter blocks.
//
// If the parameter_block_ordering contains more than one
// elimination group and support for constrained fill-reducing
// ordering is available in the sparse linear algebra library
// (SuiteSparse version >= 4.2.0) then the fill reducing
// ordering will take it into account, otherwise it will be ignored.
static bool ReorderProgramForSparseNormalCholesky(
const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
const ParameterBlockOrdering* parameter_block_ordering,
Program* program,
string* message);
// Schur type solvers require that all parameter blocks eliminated
// by the Schur eliminator occur before others and the residuals be
// sorted in lexicographic order of their parameter blocks.
//
// If the parameter_block_ordering only contains one elimination
// group then a maximal independent set is computed and used as the
// first elimination group, otherwise the user's ordering is used.
//
// If the linear solver type is SPARSE_SCHUR and support for
// constrained fill-reducing ordering is available in the sparse
// linear algebra library (SuiteSparse version >= 4.2.0) then
// columns of the schur complement matrix are ordered to reduce the
// fill-in the Cholesky factorization.
//
// Upon return, ordering contains the parameter block ordering that
// was used to order the program.
static bool ReorderProgramForSchurTypeLinearSolver(
const LinearSolverType linear_solver_type,
const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
const ProblemImpl::ParameterMap& parameter_map,
ParameterBlockOrdering* parameter_block_ordering,
Program* program,
string* message);
// array contains a list of (possibly repeating) non-negative
// integers. Let us assume that we have constructed another array
// `p` by sorting and uniqueing the entries of array.
// CompactifyArray replaces each entry in "array" with its position
// in `p`.
static void CompactifyArray(vector<int>* array);
};
} // namespace internal
} // namespace ceres
#endif // CERES_INTERNAL_SOLVER_IMPL_H_

78
extern/libmv/third_party/ceres/internal/ceres/solver_utils.cc vendored Normal file
View File

@ -0,0 +1,78 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#include <string>
#include "ceres/internal/port.h"
#include "ceres/version.h"
namespace ceres {
namespace internal {
string VersionString() {
string value = string(CERES_VERSION_STRING);
#ifdef CERES_NO_LAPACK
value += "-no_lapack";
#else
value += "-lapack";
#endif
#ifndef CERES_NO_SUITESPARSE
value += "-suitesparse";
#endif
#ifndef CERES_NO_CXSPARSE
value += "-cxsparse";
#endif
#ifdef CERES_USE_EIGEN_SPARSE
value += "-eigensparse";
#endif
#ifdef CERES_RESTRUCT_SCHUR_SPECIALIZATIONS
value += "-no_schur_specializations";
#endif
#ifdef CERES_USE_OPENMP
value += "-openmp";
#else
value += "-no_openmp";
#endif
#ifdef CERES_NO_CUSTOM_BLAS
value += "-no_custom_blas";
#endif
return value;
}
} // namespace internal
} // namespace ceres

58
extern/libmv/third_party/ceres/internal/ceres/solver_utils.h vendored Normal file
View File

@ -0,0 +1,58 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#include <algorithm>
#include <string>
namespace ceres {
namespace internal {
template <typename SummaryType>
bool IsSolutionUsable(const SummaryType& summary) {
return (summary.termination_type == CONVERGENCE ||
summary.termination_type == NO_CONVERGENCE ||
summary.termination_type == USER_SUCCESS);
}
template <typename SummaryType>
void SetSummaryFinalCost(SummaryType* summary) {
summary->final_cost = summary->initial_cost;
// We need the loop here, instead of just looking at the last
// iteration because the minimizer maybe making non-monotonic steps.
for (int i = 0; i < summary->iterations.size(); ++i) {
const IterationSummary& iteration_summary = summary->iterations[i];
summary->final_cost = min(iteration_summary.cost, summary->final_cost);
}
}
string VersionString();
} // namespace internal
} // namespace ceres

250
extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.cc vendored
View File

@ -28,11 +28,6 @@
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
// This include must come before any #ifndef check on Ceres compile options.
#include "ceres/internal/port.h"
#if !defined(CERES_NO_SUITESPARSE) || !defined(CERES_NO_CXSPARSE)
#include "ceres/sparse_normal_cholesky_solver.h"
#include <algorithm>
@ -48,9 +43,67 @@
#include "ceres/triplet_sparse_matrix.h"
#include "ceres/types.h"
#include "ceres/wall_time.h"
#include "Eigen/SparseCore"
#ifdef CERES_USE_EIGEN_SPARSE
#include "Eigen/SparseCholesky"
#endif
namespace ceres {
namespace internal {
namespace {
#ifdef CERES_USE_EIGEN_SPARSE
// A templated factorized and solve function, which allows us to use
// the same code independent of whether a AMD or a Natural ordering is
// used.
template <typename SimplicialCholeskySolver>
LinearSolver::Summary SimplicialLDLTSolve(
Eigen::MappedSparseMatrix<double, Eigen::ColMajor>& lhs,
const bool do_symbolic_analysis,
SimplicialCholeskySolver* solver,
double* rhs_and_solution,
EventLogger* event_logger) {
LinearSolver::Summary summary;
summary.num_iterations = 1;
summary.termination_type = LINEAR_SOLVER_SUCCESS;
summary.message = "Success.";
if (do_symbolic_analysis) {
solver->analyzePattern(lhs);
event_logger->AddEvent("Analyze");
if (solver->info() != Eigen::Success) {
summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
summary.message =
"Eigen failure. Unable to find symbolic factorization.";
return summary;
}
}
solver->factorize(lhs);
event_logger->AddEvent("Factorize");
if (solver->info() != Eigen::Success) {
summary.termination_type = LINEAR_SOLVER_FAILURE;
summary.message = "Eigen failure. Unable to find numeric factorization.";
return summary;
}
const Vector rhs = VectorRef(rhs_and_solution, lhs.cols());
VectorRef(rhs_and_solution, lhs.cols()) = solver->solve(rhs);
event_logger->AddEvent("Solve");
if (solver->info() != Eigen::Success) {
summary.termination_type = LINEAR_SOLVER_FAILURE;
summary.message = "Eigen failure. Unable to do triangular solve.";
return summary;
}
return summary;
}
#endif // CERES_USE_EIGEN_SPARSE
} // namespace
SparseNormalCholeskySolver::SparseNormalCholeskySolver(
const LinearSolver::Options& options)
@ -60,19 +113,15 @@ SparseNormalCholeskySolver::SparseNormalCholeskySolver(
}
void SparseNormalCholeskySolver::FreeFactorization() {
#ifndef CERES_NO_SUITESPARSE
if (factor_ != NULL) {
ss_.Free(factor_);
factor_ = NULL;
}
#endif // CERES_NO_SUITESPARSE
#ifndef CERES_NO_CXSPARSE
if (cxsparse_factor_ != NULL) {
cxsparse_.Free(cxsparse_factor_);
cxsparse_factor_ = NULL;
}
#endif // CERES_NO_CXSPARSE
}
SparseNormalCholeskySolver::~SparseNormalCholeskySolver() {
@ -111,6 +160,9 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImpl(
case CX_SPARSE:
summary = SolveImplUsingCXSparse(A, per_solve_options, x);
break;
case EIGEN_SPARSE:
summary = SolveImplUsingEigen(A, per_solve_options, x);
break;
default:
LOG(FATAL) << "Unknown sparse linear algebra library : "
<< options_.sparse_linear_algebra_library_type;
@ -123,11 +175,120 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImpl(
return summary;
}
#ifndef CERES_NO_CXSPARSE
LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingEigen(
CompressedRowSparseMatrix* A,
const LinearSolver::PerSolveOptions& per_solve_options,
double * rhs_and_solution) {
#ifndef CERES_USE_EIGEN_SPARSE
LinearSolver::Summary summary;
summary.num_iterations = 0;
summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
summary.message =
"SPARSE_NORMAL_CHOLESKY cannot be used with EIGEN_SPARSE "
"because Ceres was not built with support for "
"Eigen's SimplicialLDLT decomposition. "
"This requires enabling building with -DEIGENSPARSE=ON.";
return summary;
#else
EventLogger event_logger("SparseNormalCholeskySolver::Eigen::Solve");
// Compute the normal equations. J'J delta = J'f and solve them
// using a sparse Cholesky factorization. Notice that when compared
// to SuiteSparse we have to explicitly compute the normal equations
// before they can be factorized. CHOLMOD/SuiteSparse on the other
// hand can just work off of Jt to compute the Cholesky
// factorization of the normal equations.
//
// TODO(sameeragarwal): See note about how this maybe a bad idea for
// dynamic sparsity.
if (outer_product_.get() == NULL) {
outer_product_.reset(
CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
*A, &pattern_));
}
CompressedRowSparseMatrix::ComputeOuterProduct(
*A, pattern_, outer_product_.get());
// Map to an upper triangular column major matrix.
//
// outer_product_ is a compressed row sparse matrix and in lower
// triangular form, when mapped to a compressed column sparse
// matrix, it becomes an upper triangular matrix.
Eigen::MappedSparseMatrix<double, Eigen::ColMajor> AtA(
outer_product_->num_rows(),
outer_product_->num_rows(),
outer_product_->num_nonzeros(),
outer_product_->mutable_rows(),
outer_product_->mutable_cols(),
outer_product_->mutable_values());
// Dynamic sparsity means that we cannot depend on a static analysis
// of sparsity structure of the jacobian, so we compute a new
// symbolic factorization every time.
if (options_.dynamic_sparsity) {
amd_ldlt_.reset(NULL);
}
bool do_symbolic_analysis = false;
// If using post ordering or dynamic sparsity, or an old version of
// Eigen, we cannot depend on a preordered jacobian, so we work with
// a SimplicialLDLT decomposition with AMD ordering.
if (options_.use_postordering ||
options_.dynamic_sparsity ||
!EIGEN_VERSION_AT_LEAST(3, 2, 2)) {
if (amd_ldlt_.get() == NULL) {
amd_ldlt_.reset(new SimplicialLDLTWithAMDOrdering);
do_symbolic_analysis = true;
}
return SimplicialLDLTSolve(AtA,
do_symbolic_analysis,
amd_ldlt_.get(),
rhs_and_solution,
&event_logger);
}
#if EIGEN_VERSION_AT_LEAST(3,2,2)
// The common case
if (natural_ldlt_.get() == NULL) {
natural_ldlt_.reset(new SimplicialLDLTWithNaturalOrdering);
do_symbolic_analysis = true;
}
return SimplicialLDLTSolve(AtA,
do_symbolic_analysis,
natural_ldlt_.get(),
rhs_and_solution,
&event_logger);
#endif
#endif // EIGEN_USE_EIGEN_SPARSE
}
LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse(
CompressedRowSparseMatrix* A,
const LinearSolver::PerSolveOptions& per_solve_options,
double * rhs_and_solution) {
#ifdef CERES_NO_CXSPARSE
LinearSolver::Summary summary;
summary.num_iterations = 0;
summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
summary.message =
"SPARSE_NORMAL_CHOLESKY cannot be used with CX_SPARSE "
"because Ceres was not built with support for CXSparse. "
"This requires enabling building with -DCXSPARSE=ON.";
return summary;
#else
EventLogger event_logger("SparseNormalCholeskySolver::CXSparse::Solve");
LinearSolver::Summary summary;
@ -137,11 +298,14 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse(
// Compute the normal equations. J'J delta = J'f and solve them
// using a sparse Cholesky factorization. Notice that when compared
// to SuiteSparse we have to explicitly compute the transpose of Jt,
// and then the normal equations before they can be
// factorized. CHOLMOD/SuiteSparse on the other hand can just work
// off of Jt to compute the Cholesky factorization of the normal
// equations.
// to SuiteSparse we have to explicitly compute the normal equations
// before they can be factorized. CHOLMOD/SuiteSparse on the other
// hand can just work off of Jt to compute the Cholesky
// factorization of the normal equations.
//
// TODO(sameeragarwal): If dynamic sparsity is enabled, then this is
// not a good idea performance wise, since the jacobian has far too
// many entries and the program will go crazy with memory.
if (outer_product_.get() == NULL) {
outer_product_.reset(
CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
@ -179,30 +343,35 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse(
summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
summary.message =
"CXSparse failure. Unable to find symbolic factorization.";
} else if (!cxsparse_.SolveCholesky(AtA, cxsparse_factor_, rhs_and_solution)) {
} else if (!cxsparse_.SolveCholesky(AtA,
cxsparse_factor_,
rhs_and_solution)) {
summary.termination_type = LINEAR_SOLVER_FAILURE;
summary.message = "CXSparse::SolveCholesky failed.";
}
event_logger.AddEvent("Solve");
return summary;
}
#else
LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse(
CompressedRowSparseMatrix* A,
const LinearSolver::PerSolveOptions& per_solve_options,
double * rhs_and_solution) {
LOG(FATAL) << "No CXSparse support in Ceres.";
// Unreachable but MSVC does not know this.
return LinearSolver::Summary();
}
#endif
}
#ifndef CERES_NO_SUITESPARSE
LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse(
CompressedRowSparseMatrix* A,
const LinearSolver::PerSolveOptions& per_solve_options,
double * rhs_and_solution) {
#ifdef CERES_NO_SUITESPARSE
LinearSolver::Summary summary;
summary.num_iterations = 0;
summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
summary.message =
"SPARSE_NORMAL_CHOLESKY cannot be used with SUITE_SPARSE "
"because Ceres was not built with support for SuiteSparse. "
"This requires enabling building with -DSUITESPARSE=ON.";
return summary;
#else
EventLogger event_logger("SparseNormalCholeskySolver::SuiteSparse::Solve");
LinearSolver::Summary summary;
summary.termination_type = LINEAR_SOLVER_SUCCESS;
@ -226,7 +395,8 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse(
if (options_.dynamic_sparsity) {
factor_ = ss_.AnalyzeCholesky(&lhs, &summary.message);
} else {
factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs, &summary.message);
factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs,
&summary.message);
}
}
}
@ -234,6 +404,8 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse(
if (factor_ == NULL) {
summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
// No need to set message as it has already been set by the
// symbolic analysis routines above.
return summary;
}
@ -242,7 +414,9 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse(
return summary;
}
cholmod_dense* rhs = ss_.CreateDenseVector(rhs_and_solution, num_cols, num_cols);
cholmod_dense* rhs = ss_.CreateDenseVector(rhs_and_solution,
num_cols,
num_cols);
cholmod_dense* solution = ss_.Solve(factor_, rhs, &summary.message);
event_logger.AddEvent("Solve");
@ -251,25 +425,15 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse(
memcpy(rhs_and_solution, solution->x, num_cols * sizeof(*rhs_and_solution));
ss_.Free(solution);
} else {
// No need to set message as it has already been set by the
// numeric factorization routine above.
summary.termination_type = LINEAR_SOLVER_FAILURE;
}
event_logger.AddEvent("Teardown");
return summary;
}
#else
LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse(
CompressedRowSparseMatrix* A,
const LinearSolver::PerSolveOptions& per_solve_options,
double * rhs_and_solution) {
LOG(FATAL) << "No SuiteSparse support in Ceres.";
// Unreachable but MSVC does not know this.
return LinearSolver::Summary();
}
#endif
}
} // namespace internal
} // namespace ceres
#endif // !defined(CERES_NO_SUITESPARSE) || !defined(CERES_NO_CXSPARSE)

43
extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.h vendored
View File

@ -34,15 +34,19 @@
#ifndef CERES_INTERNAL_SPARSE_NORMAL_CHOLESKY_SOLVER_H_
#define CERES_INTERNAL_SPARSE_NORMAL_CHOLESKY_SOLVER_H_
#include <vector>
// This include must come before any #ifndef check on Ceres compile options.
#include "ceres/internal/port.h"
#if !defined(CERES_NO_SUITESPARSE) || !defined(CERES_NO_CXSPARSE)
#include "ceres/cxsparse.h"
#include "ceres/internal/macros.h"
#include "ceres/linear_solver.h"
#include "ceres/suitesparse.h"
#include "ceres/cxsparse.h"
#ifdef CERES_USE_EIGEN_SPARSE
#include "Eigen/SparseCholesky"
#endif
namespace ceres {
namespace internal {
@ -74,6 +78,12 @@ class SparseNormalCholeskySolver : public CompressedRowSparseMatrixSolver {
const LinearSolver::PerSolveOptions& options,
double* rhs_and_solution);
// Crashes if CERES_USE_EIGEN_SPARSE is not defined.
LinearSolver::Summary SolveImplUsingEigen(
CompressedRowSparseMatrix* A,
const LinearSolver::PerSolveOptions& options,
double* rhs_and_solution);
void FreeFactorization();
SuiteSparse ss_;
@ -83,6 +93,32 @@ class SparseNormalCholeskySolver : public CompressedRowSparseMatrixSolver {
CXSparse cxsparse_;
// Cached factorization
cs_dis* cxsparse_factor_;
#ifdef CERES_USE_EIGEN_SPARSE
// The preprocessor gymnastics here are dealing with the fact that
// before version 3.2.2, Eigen did not support a third template
// parameter to specify the ordering.
#if EIGEN_VERSION_AT_LEAST(3,2,2)
typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>, Eigen::Upper,
Eigen::NaturalOrdering<int> >
SimplicialLDLTWithNaturalOrdering;
scoped_ptr<SimplicialLDLTWithNaturalOrdering> natural_ldlt_;
typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>, Eigen::Upper,
Eigen::AMDOrdering<int> >
SimplicialLDLTWithAMDOrdering;
scoped_ptr<SimplicialLDLTWithAMDOrdering> amd_ldlt_;
#else
typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>, Eigen::Upper>
SimplicialLDLTWithAMDOrdering;
scoped_ptr<SimplicialLDLTWithAMDOrdering> amd_ldlt_;
#endif
#endif
scoped_ptr<CompressedRowSparseMatrix> outer_product_;
vector<int> pattern_;
const LinearSolver::Options options_;
@ -92,5 +128,4 @@ class SparseNormalCholeskySolver : public CompressedRowSparseMatrixSolver {
} // namespace internal
} // namespace ceres
#endif // !defined(CERES_NO_SUITESPARSE) || !defined(CERES_NO_CXSPARSE)
#endif // CERES_INTERNAL_SPARSE_NORMAL_CHOLESKY_SOLVER_H_

17
extern/libmv/third_party/ceres/internal/ceres/suitesparse.h vendored
View File

@ -283,9 +283,24 @@ class SuiteSparse {
#else // CERES_NO_SUITESPARSE
class SuiteSparse {};
typedef void cholmod_factor;
class SuiteSparse {
public:
// Defining this static function even when SuiteSparse is not
// available, allows client code to check for the presence of CAMD
// without checking for the absence of the CERES_NO_CAMD symbol.
//
// This is safer because the symbol maybe missing due to a user
// accidently not including suitesparse.h in their code when
// checking for the symbol.
static bool IsConstrainedApproximateMinimumDegreeOrderingAvailable() {
return false;
}
void Free(void*) {};
};
#endif // CERES_NO_SUITESPARSE
#endif // CERES_INTERNAL_SUITESPARSE_H_

13
extern/libmv/third_party/ceres/internal/ceres/trust_region_minimizer.cc vendored
View File

@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2012 Google Inc. All rights reserved.
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
@ -40,6 +40,7 @@
#include "Eigen/Core"
#include "ceres/array_utils.h"
#include "ceres/coordinate_descent_minimizer.h"
#include "ceres/evaluator.h"
#include "ceres/file.h"
#include "ceres/internal/eigen.h"
@ -128,9 +129,10 @@ void TrustRegionMinimizer::Minimize(const Minimizer::Options& options,
double iteration_start_time = start_time;
Init(options);
Evaluator* evaluator = CHECK_NOTNULL(options_.evaluator);
SparseMatrix* jacobian = CHECK_NOTNULL(options_.jacobian);
TrustRegionStrategy* strategy = CHECK_NOTNULL(options_.trust_region_strategy);
Evaluator* evaluator = CHECK_NOTNULL(options_.evaluator.get());
SparseMatrix* jacobian = CHECK_NOTNULL(options_.jacobian.get());
TrustRegionStrategy* strategy =
CHECK_NOTNULL(options_.trust_region_strategy.get());
const bool is_not_silent = !options.is_silent;
@ -253,7 +255,8 @@ void TrustRegionMinimizer::Minimize(const Minimizer::Options& options,
double candidate_cost = cost;
double accumulated_candidate_model_cost_change = 0.0;
int num_consecutive_invalid_steps = 0;
bool inner_iterations_are_enabled = options.inner_iteration_minimizer != NULL;
bool inner_iterations_are_enabled =
options.inner_iteration_minimizer.get() != NULL;
while (true) {
bool inner_iterations_were_useful = false;
if (!RunCallbacks(options, iteration_summary, summary)) {

356
extern/libmv/third_party/ceres/internal/ceres/trust_region_preprocessor.cc vendored Normal file
View File

@ -0,0 +1,356 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 Google Inc. All rights reserved.
// http://code.google.com/p/ceres-solver/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#include "ceres/trust_region_preprocessor.h"
#include <numeric>
#include <string>
#include "ceres/callbacks.h"
#include "ceres/evaluator.h"
#include "ceres/linear_solver.h"
#include "ceres/minimizer.h"
#include "ceres/parameter_block.h"
#include "ceres/preconditioner.h"
#include "ceres/preprocessor.h"
#include "ceres/problem_impl.h"
#include "ceres/program.h"
#include "ceres/reorder_program.h"
#include "ceres/suitesparse.h"
#include "ceres/trust_region_strategy.h"
#include "ceres/wall_time.h"
namespace ceres {
namespace internal {
namespace {
ParameterBlockOrdering* CreateDefaultLinearSolverOrdering(
const Program& program) {
ParameterBlockOrdering* ordering = new ParameterBlockOrdering;
const vector<ParameterBlock*>& parameter_blocks = program.parameter_blocks();
for (int i = 0; i < parameter_blocks.size(); ++i) {
ordering->AddElementToGroup(
const_cast<double*>(parameter_blocks[i]->user_state()), 0);
}
return ordering;
}
// Check if all the user supplied values in the parameter blocks are
// sane or not, and if the program is feasible or not.
bool IsProgramValid(const Program& program, string* error) {
return (program.ParameterBlocksAreFinite(error) &&
program.IsFeasible(error));
}
void AlternateLinearSolverAndPreconditionerForSchurTypeLinearSolver(
Solver::Options* options) {
if (!IsSchurType(options->linear_solver_type)) {
return;
}
const LinearSolverType linear_solver_type_given = options->linear_solver_type;
const PreconditionerType preconditioner_type_given =
options->preconditioner_type;
options->linear_solver_type = LinearSolver::LinearSolverForZeroEBlocks(
linear_solver_type_given);
string message;
if (linear_solver_type_given == ITERATIVE_SCHUR) {
options->preconditioner_type = Preconditioner::PreconditionerForZeroEBlocks(
preconditioner_type_given);
message =
StringPrintf(
"No E blocks. Switching from %s(%s) to %s(%s).",
LinearSolverTypeToString(linear_solver_type_given),
PreconditionerTypeToString(preconditioner_type_given),
LinearSolverTypeToString(options->linear_solver_type),
PreconditionerTypeToString(options->preconditioner_type));
} else {
message =
StringPrintf(
"No E blocks. Switching from %s to %s.",
LinearSolverTypeToString(linear_solver_type_given),
LinearSolverTypeToString(options->linear_solver_type));
}
VLOG_IF(1, options->logging_type != SILENT) << message;
}
// For Schur type and SPARSE_NORMAL_CHOLESKY linear solvers, reorder
// the program to reduce fill-in and increase cache coherency.
bool ReorderProgram(PreprocessedProblem* pp) {
Solver::Options& options = pp->options;
if (IsSchurType(options.linear_solver_type)) {
return ReorderProgramForSchurTypeLinearSolver(
options.linear_solver_type,
options.sparse_linear_algebra_library_type,
pp->problem->parameter_map(),
options.linear_solver_ordering.get(),
pp->reduced_program.get(),
&pp->error);
}
if (options.linear_solver_type == SPARSE_NORMAL_CHOLESKY &&
!options.dynamic_sparsity) {
return ReorderProgramForSparseNormalCholesky(
options.sparse_linear_algebra_library_type,
*options.linear_solver_ordering,
pp->reduced_program.get(),
&pp->error);
}
return true;
}
// Configure and create a linear solver object. In doing so, if a
// sparse direct factorization based linear solver is being used, then
// find a fill reducing ordering and reorder the program as needed
// too.
bool SetupLinearSolver(PreprocessedProblem* pp) {
Solver::Options& options = pp->options;
if (options.linear_solver_ordering.get() == NULL) {
// If the user has not supplied a linear solver ordering, then we
// assume that they are giving all the freedom to us in choosing
// the best possible ordering. This intent can be indicated by
// putting all the parameter blocks in the same elimination group.
options.linear_solver_ordering.reset(
CreateDefaultLinearSolverOrdering(*pp->reduced_program));
} else {
// If the user supplied an ordering, then check if the first
// elimination group is still non-empty after the reduced problem
// has been constructed.
//
// This is important for Schur type linear solvers, where the
// first elimination group is special -- it needs to be an
// independent set.
//
// If the first elimination group is empty, then we cannot use the
// user's requested linear solver (and a preconditioner as the
// case may be) so we must use a different one.
ParameterBlockOrdering* ordering = options.linear_solver_ordering.get();
const int min_group_id = ordering->MinNonZeroGroup();
ordering->Remove(pp->removed_parameter_blocks);
if (IsSchurType(options.linear_solver_type) &&
min_group_id != ordering->MinNonZeroGroup()) {
AlternateLinearSolverAndPreconditionerForSchurTypeLinearSolver(
&options);
}
}
// Reorder the program to reduce fill in and improve cache coherency
// of the Jacobian.
if (!ReorderProgram(pp)) {
return false;
}
// Configure the linear solver.
pp->linear_solver_options = LinearSolver::Options();
pp->linear_solver_options.min_num_iterations =
options.min_linear_solver_iterations;
pp->linear_solver_options.max_num_iterations =
options.max_linear_solver_iterations;
pp->linear_solver_options.type = options.linear_solver_type;
pp->linear_solver_options.preconditioner_type = options.preconditioner_type;
pp->linear_solver_options.visibility_clustering_type =
options.visibility_clustering_type;
pp->linear_solver_options.sparse_linear_algebra_library_type =
options.sparse_linear_algebra_library_type;
pp->linear_solver_options.dense_linear_algebra_library_type =
options.dense_linear_algebra_library_type;
pp->linear_solver_options.dynamic_sparsity = options.dynamic_sparsity;
pp->linear_solver_options.num_threads = options.num_linear_solver_threads;
// Ignore user's postordering preferences and force it to be true if
// cholmod_camd is not available. This ensures that the linear
// solver does not assume that a fill-reducing pre-ordering has been
// done.
pp->linear_solver_options.use_postordering = options.use_postordering;
if (options.linear_solver_type == SPARSE_SCHUR &&
options.sparse_linear_algebra_library_type == SUITE_SPARSE &&
!SuiteSparse::IsConstrainedApproximateMinimumDegreeOrderingAvailable()) {
pp->linear_solver_options.use_postordering = true;
}
OrderingToGroupSizes(options.linear_solver_ordering.get(),
&pp->linear_solver_options.elimination_groups);
// Schur type solvers expect at least two elimination groups. If
// there is only one elimination group, then it is guaranteed that
// this group only contains e_blocks. Thus we add a dummy
// elimination group with zero blocks in it.
if (IsSchurType(pp->linear_solver_options.type) &&
pp->linear_solver_options.elimination_groups.size() == 1) {
pp->linear_solver_options.elimination_groups.push_back(0);
}
pp->linear_solver.reset(LinearSolver::Create(pp->linear_solver_options));
return (pp->linear_solver.get() != NULL);
}
// Configure and create the evaluator.
bool SetupEvaluator(PreprocessedProblem* pp) {
const Solver::Options& options = pp->options;
pp->evaluator_options = Evaluator::Options();
pp->evaluator_options.linear_solver_type = options.linear_solver_type;
pp->evaluator_options.num_eliminate_blocks = 0;
if (IsSchurType(options.linear_solver_type)) {
pp->evaluator_options.num_eliminate_blocks =
options
.linear_solver_ordering
->group_to_elements().begin()
->second.size();
}
pp->evaluator_options.num_threads = options.num_threads;
pp->evaluator_options.dynamic_sparsity = options.dynamic_sparsity;
pp->evaluator.reset(Evaluator::Create(pp->evaluator_options,
pp->reduced_program.get(),
&pp->error));
return (pp->evaluator.get() != NULL);
}
// If the user requested inner iterations, then find an inner
// iteration ordering as needed and configure and create a
// CoordinateDescentMinimizer object to perform the inner iterations.
bool SetupInnerIterationMinimizer(PreprocessedProblem* pp) {
Solver::Options& options = pp->options;
if (!options.use_inner_iterations) {
return true;
}
// With just one parameter block, the outer iteration of the trust
// region method and inner iterations are doing exactly the same
// thing, and thus inner iterations are not needed.
if (pp->reduced_program->NumParameterBlocks() == 1) {
LOG(WARNING) << "Reduced problem only contains one parameter block."
<< "Disabling inner iterations.";
return true;
}
if (options.inner_iteration_ordering.get() != NULL) {
// If the user supplied an ordering, then remove the set of
// inactive parameter blocks from it
options.inner_iteration_ordering->Remove(pp->removed_parameter_blocks);
if (options.inner_iteration_ordering->NumElements() == 0) {
LOG(WARNING) << "No remaining elements in the inner iteration ordering.";
return true;
}
// Validate the reduced ordering.
if (!CoordinateDescentMinimizer::IsOrderingValid(
*pp->reduced_program,
*options.inner_iteration_ordering,
&pp->error)) {
return false;
}
} else {
// The user did not supply an ordering, so create one.
options.inner_iteration_ordering.reset(
CoordinateDescentMinimizer::CreateOrdering(*pp->reduced_program));
}
pp->inner_iteration_minimizer.reset(new CoordinateDescentMinimizer);
return pp->inner_iteration_minimizer->Init(*pp->reduced_program,
pp->problem->parameter_map(),
*options.inner_iteration_ordering,
&pp->error);
}
// Configure and create a TrustRegionMinimizer object.
void SetupMinimizerOptions(PreprocessedProblem* pp) {
const Solver::Options& options = pp->options;
SetupCommonMinimizerOptions(pp);
pp->minimizer_options.is_constrained =
pp->reduced_program->IsBoundsConstrained();
pp->minimizer_options.jacobian.reset(pp->evaluator->CreateJacobian());
pp->minimizer_options.inner_iteration_minimizer =
pp->inner_iteration_minimizer;
TrustRegionStrategy::Options strategy_options;
strategy_options.linear_solver = pp->linear_solver.get();
strategy_options.initial_radius =
options.initial_trust_region_radius;
strategy_options.max_radius = options.max_trust_region_radius;
strategy_options.min_lm_diagonal = options.min_lm_diagonal;
strategy_options.max_lm_diagonal = options.max_lm_diagonal;
strategy_options.trust_region_strategy_type =
options.trust_region_strategy_type;
strategy_options.dogleg_type = options.dogleg_type;
pp->minimizer_options.trust_region_strategy.reset(
CHECK_NOTNULL(TrustRegionStrategy::Create(strategy_options)));
}
} // namespace
TrustRegionPreprocessor::~TrustRegionPreprocessor() {
}
bool TrustRegionPreprocessor::Preprocess(const Solver::Options& options,
ProblemImpl* problem,
PreprocessedProblem* pp) {
CHECK_NOTNULL(pp);
pp->options = options;
ChangeNumThreadsIfNeeded(&pp->options);
pp->problem = problem;
Program* program = problem->mutable_program();
if (!IsProgramValid(*program, &pp->error)) {
return false;
}
pp->reduced_program.reset(
program->CreateReducedProgram(&pp->removed_parameter_blocks,
&pp->fixed_cost,
&pp->error));
if (pp->reduced_program.get() == NULL) {
return false;
}
if (pp->reduced_program->NumParameterBlocks() == 0) {
// The reduced problem has no parameter or residual blocks. There
// is nothing more to do.
return true;
}
if (!SetupLinearSolver(pp) ||
!SetupEvaluator(pp) ||
!SetupInnerIterationMinimizer(pp)) {
return false;
}
SetupMinimizerOptions(pp);
return true;
}
} // namespace internal
} // namespace ceres

50
extern/libmv/third_party/ceres/internal/ceres/trust_region_preprocessor.h vendored Normal file
View File

@ -0,0 +1,50 @@
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2014 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: sameragarwal@google.com (Sameer Agarwal)
#ifndef CERES_INTERNAL_TRUST_REGION_PREPROCESSOR_H_
#define CERES_INTERNAL_TRUST_REGION_PREPROCESSOR_H_
#include "ceres/preprocessor.h"
namespace ceres {
namespace internal {
class TrustRegionPreprocessor : public Preprocessor {
public:
virtual ~TrustRegionPreprocessor();
virtual bool Preprocess(const Solver::Options& options,
ProblemImpl* problem,
PreprocessedProblem* preprocessed_problem);
};
} // namespace internal
} // namespace ceres
#endif // CERES_INTERNAL_TRUST_REGION_PREPROCESSOR_H_

14
extern/libmv/third_party/ceres/internal/ceres/types.cc vendored
View File

@ -96,6 +96,7 @@ const char* SparseLinearAlgebraLibraryTypeToString(
switch (type) {
CASESTR(SUITE_SPARSE);
CASESTR(CX_SPARSE);
CASESTR(EIGEN_SPARSE);
default:
return "UNKNOWN";
}
@ -107,6 +108,7 @@ bool StringToSparseLinearAlgebraLibraryType(
UpperCase(&value);
STRENUM(SUITE_SPARSE);
STRENUM(CX_SPARSE);
STRENUM(EIGEN_SPARSE);
return false;
}
@ -240,7 +242,7 @@ const char* NonlinearConjugateGradientTypeToString(
NonlinearConjugateGradientType type) {
switch (type) {
CASESTR(FLETCHER_REEVES);
CASESTR(POLAK_RIBIRERE);
CASESTR(POLAK_RIBIERE);
CASESTR(HESTENES_STIEFEL);
default:
return "UNKNOWN";
@ -252,7 +254,7 @@ bool StringToNonlinearConjugateGradientType(
NonlinearConjugateGradientType* type) {
UpperCase(&value);
STRENUM(FLETCHER_REEVES);
STRENUM(POLAK_RIBIRERE);
STRENUM(POLAK_RIBIERE);
STRENUM(HESTENES_STIEFEL);
return false;
}
@ -261,8 +263,8 @@ const char* CovarianceAlgorithmTypeToString(
CovarianceAlgorithmType type) {
switch (type) {
CASESTR(DENSE_SVD);
CASESTR(SPARSE_CHOLESKY);
CASESTR(SPARSE_QR);
CASESTR(EIGEN_SPARSE_QR);
CASESTR(SUITE_SPARSE_QR);
default:
return "UNKNOWN";
}
@ -273,8 +275,8 @@ bool StringToCovarianceAlgorithmType(
CovarianceAlgorithmType* type) {
UpperCase(&value);
STRENUM(DENSE_SVD);
STRENUM(SPARSE_CHOLESKY);
STRENUM(SPARSE_QR);
STRENUM(EIGEN_SPARSE_QR);
STRENUM(SUITE_SPARSE_QR);
return false;
}

5
extern/libmv/third_party/ceres/internal/ceres/visibility.cc vendored
View File

@ -76,7 +76,8 @@ void ComputeVisibility(const CompressedRowBlockStructure& block_structure,
}
}
Graph<int>* CreateSchurComplementGraph(const vector<set<int> >& visibility) {
WeightedGraph<int>* CreateSchurComplementGraph(
const vector<set<int> >& visibility) {
const time_t start_time = time(NULL);
// Compute the number of e_blocks/point blocks. Since the visibility
// set for each e_block/camera contains the set of e_blocks/points
@ -122,7 +123,7 @@ Graph<int>* CreateSchurComplementGraph(const vector<set<int> >& visibility) {
}
}
Graph<int>* graph = new Graph<int>();
WeightedGraph<int>* graph = new WeightedGraph<int>;
// Add vertices and initialize the pairs for self edges so that self
// edges are guaranteed. This is needed for the Canonical views

5
extern/libmv/third_party/ceres/internal/ceres/visibility.h vendored
View File

@ -72,9 +72,10 @@ void ComputeVisibility(const CompressedRowBlockStructure& block_structure,
// matrix/Schur complement matrix obtained by eliminating the e_blocks
// from the normal equations.
//
// Caller acquires ownership of the returned Graph pointer
// Caller acquires ownership of the returned WeightedGraph pointer
// (heap-allocated).
Graph<int>* CreateSchurComplementGraph(const vector<set<int> >& visibility);
WeightedGraph<int>* CreateSchurComplementGraph(
const vector<set<int> >& visibility);
} // namespace internal
} // namespace ceres

12
extern/libmv/third_party/ceres/internal/ceres/visibility_based_preconditioner.cc vendored
View File

@ -167,9 +167,9 @@ void VisibilityBasedPreconditioner::ComputeClusterTridiagonalSparsity(
// maximum spanning forest of this graph.
vector<set<int> > cluster_visibility;
ComputeClusterVisibility(visibility, &cluster_visibility);
scoped_ptr<Graph<int> > cluster_graph(
scoped_ptr<WeightedGraph<int> > cluster_graph(
CHECK_NOTNULL(CreateClusterGraph(cluster_visibility)));
scoped_ptr<Graph<int> > forest(
scoped_ptr<WeightedGraph<int> > forest(
CHECK_NOTNULL(Degree2MaximumSpanningForest(*cluster_graph)));
ForestToClusterPairs(*forest, &cluster_pairs_);
}
@ -189,7 +189,7 @@ void VisibilityBasedPreconditioner::InitStorage(
// memberships for each camera block.
void VisibilityBasedPreconditioner::ClusterCameras(
const vector<set<int> >& visibility) {
scoped_ptr<Graph<int> > schur_complement_graph(
scoped_ptr<WeightedGraph<int> > schur_complement_graph(
CHECK_NOTNULL(CreateSchurComplementGraph(visibility)));
HashMap<int, int> membership;
@ -498,7 +498,7 @@ bool VisibilityBasedPreconditioner::IsBlockPairOffDiagonal(
// Convert a graph into a list of edges that includes self edges for
// each vertex.
void VisibilityBasedPreconditioner::ForestToClusterPairs(
const Graph<int>& forest,
const WeightedGraph<int>& forest,
HashSet<pair<int, int> >* cluster_pairs) const {
CHECK_NOTNULL(cluster_pairs)->clear();
const HashSet<int>& vertices = forest.vertices();
@ -541,9 +541,9 @@ void VisibilityBasedPreconditioner::ComputeClusterVisibility(
// Construct a graph whose vertices are the clusters, and the edge
// weights are the number of 3D points visible to cameras in both the
// vertices.
Graph<int>* VisibilityBasedPreconditioner::CreateClusterGraph(
WeightedGraph<int>* VisibilityBasedPreconditioner::CreateClusterGraph(
const vector<set<int> >& cluster_visibility) const {
Graph<int>* cluster_graph = new Graph<int>;
WeightedGraph<int>* cluster_graph = new WeightedGraph<int>;
for (int i = 0; i < num_clusters_; ++i) {
cluster_graph->AddVertex(i);

5
extern/libmv/third_party/ceres/internal/ceres/visibility_based_preconditioner.h vendored
View File

@ -156,8 +156,9 @@ class VisibilityBasedPreconditioner : public BlockSparseMatrixPreconditioner {
vector<int>* membership_vector) const;
void ComputeClusterVisibility(const vector<set<int> >& visibility,
vector<set<int> >* cluster_visibility) const;
Graph<int>* CreateClusterGraph(const vector<set<int> >& visibility) const;
void ForestToClusterPairs(const Graph<int>& forest,
WeightedGraph<int>* CreateClusterGraph(
const vector<set<int> >& visibility) const;
void ForestToClusterPairs(const WeightedGraph<int>& forest,
HashSet<pair<int, int> >* cluster_pairs) const;
void ComputeBlockPairsInPreconditioner(const CompressedRowBlockStructure& bs);
bool IsBlockPairInPreconditioner(int block1, int block2) const;