Tomato: re-integrate ceres library with updates needed for tracking improvement
authorSergey Sharybin <sergey.vfx@gmail.com>
Thu, 10 May 2012 10:17:59 +0000 (10:17 +0000)
committerSergey Sharybin <sergey.vfx@gmail.com>
Thu, 10 May 2012 10:17:59 +0000 (10:17 +0000)
Also commit missed patch.

21 files changed:
extern/libmv/third_party/ceres/CMakeLists.txt
extern/libmv/third_party/ceres/ChangeLog
extern/libmv/third_party/ceres/include/ceres/autodiff_cost_function.h
extern/libmv/third_party/ceres/include/ceres/internal/autodiff.h
extern/libmv/third_party/ceres/include/ceres/sized_cost_function.h
extern/libmv/third_party/ceres/include/ceres/solver.h
extern/libmv/third_party/ceres/include/ceres/types.h
extern/libmv/third_party/ceres/internal/ceres/block_sparse_matrix.cc
extern/libmv/third_party/ceres/internal/ceres/block_sparse_matrix.h
extern/libmv/third_party/ceres/internal/ceres/compressed_row_sparse_matrix.cc
extern/libmv/third_party/ceres/internal/ceres/compressed_row_sparse_matrix.h
extern/libmv/third_party/ceres/internal/ceres/dense_sparse_matrix.cc
extern/libmv/third_party/ceres/internal/ceres/dense_sparse_matrix.h
extern/libmv/third_party/ceres/internal/ceres/levenberg_marquardt.cc
extern/libmv/third_party/ceres/internal/ceres/linear_least_squares_problems.cc
extern/libmv/third_party/ceres/internal/ceres/linear_least_squares_problems.h
extern/libmv/third_party/ceres/internal/ceres/minimizer.h
extern/libmv/third_party/ceres/internal/ceres/sparse_matrix.h
extern/libmv/third_party/ceres/internal/ceres/triplet_sparse_matrix.cc
extern/libmv/third_party/ceres/internal/ceres/triplet_sparse_matrix.h
extern/libmv/third_party/ceres/patches/msvc_isfinite.patch [new file with mode: 0644]

index 19158dcfebb07b2471f82d080f10bc8f1e9c2aa7..5207bddec12558f3e0d4e368661365556db86186 100644 (file)
@@ -208,11 +208,11 @@ else()
 endif()
 
 add_definitions(
-    -DCERES_HAVE_PTHREAD
-    -D"CERES_HASH_NAMESPACE_START=namespace std { namespace tr1 {"
-    -D"CERES_HASH_NAMESPACE_END=}}"
-    -DCERES_NO_SUITESPARSE
-    -DCERES_DONT_HAVE_PROTOCOL_BUFFERS
+       -DCERES_HAVE_PTHREAD
+       -D"CERES_HASH_NAMESPACE_START=namespace std { namespace tr1 {"
+       -D"CERES_HASH_NAMESPACE_END=}}"
+       -DCERES_NO_SUITESPARSE
+       -DCERES_DONT_HAVE_PROTOCOL_BUFFERS
 )
 
 blender_add_lib(extern_ceres "${SRC}" "${INC}" "${INC_SYS}")
index c652a73c0c180b38326c924fdfd5e8a1d2af5188..6e919658f13598612d561a5c5f3675e076c95682 100644 (file)
@@ -1,3 +1,100 @@
+commit ca72152362ae1f4b9928c012e74b4d49d094a4ca
+Merge: d297f8d 0a04199
+Author: Keir Mierle <mierle@gmail.com>
+Date:   Wed May 9 13:10:59 2012 -0700
+
+    Merge branch 'master' into windows
+
+commit 0a04199ef279cc9ea97f665fed8e7fae717813c3
+Merge: fdeb577 f2571f1
+Author: Keir Mierle <mierle@gmail.com>
+Date:   Wed May 9 12:54:56 2012 -0700
+
+    Merge branch 'master' of https://code.google.com/p/ceres-solver
+
+commit fdeb5772cc5eeebca4d776d220d80cc91b6d0f74
+Author: Keir Mierle <mierle@gmail.com>
+Date:   Wed May 9 07:38:07 2012 -0700
+
+    Support varying numbers of residuals in autodiff.
+    
+    This commit modifies the only function in autodiff that takes a
+    templated number of outputs (i.e. residuals) and makes that
+    template parameter a normal parameter. With that change, it
+    is a trivial matter to support a dynamic number of residuals.
+    
+    The API for dynamic residuals is to pass a fake number of
+    residuals as the second template argument to
+    AutoDiffCostFunction, and to pass the real number of
+    parameters as a second constructor argument.
+
+commit da3e0563cc12e08e7b3e0fbf11d9cc8cfe9658aa
+Author: Sameer Agarwal <sameeragarwal@google.com>
+Date:   Wed May 9 11:57:47 2012 -0700
+
+    Typo corrections in the documentation from Bing
+
+commit aa9526d8e8fb34c23d63e3af5bf9239b0c4ea603
+Author: Sameer Agarwal <sameeragarwal@google.com>
+Date:   Tue May 8 21:22:09 2012 -0700
+
+    Share search paths across various library searches.
+    Fix typos in glog search.
+    Split the error messages for include and lib.
+    Enable building of tests by default.
+    Made building on homebrew installations a bit better.
+    Remove temporary variables for glog and gflags.
+
+commit f2571f186850ed3dd316236ac4be488979df7d30
+Author: Sameer Agarwal <sameeragarwal@google.com>
+Date:   Wed May 9 11:57:47 2012 -0700
+
+    Typo corrections in the documentation from Bing
+
+commit 8f7f11ff7d07737435428a2620c52419cf99f98e
+Merge: e6c17c4 eaccbb3
+Author: Sameer Agarwal <sameeragarwal@google.com>
+Date:   Wed May 9 11:34:15 2012 -0700
+
+    Merge branch 'master' of https://code.google.com/p/ceres-solver
+
+commit e6c17c4c9d9307218f6f739cea39bc2d87733d4d
+Author: Sameer Agarwal <sameeragarwal@google.com>
+Date:   Tue May 8 21:22:09 2012 -0700
+
+    Share search paths across various library searches.
+    Fix typos in glog search.
+    Split the error messages for include and lib.
+    Enable building of tests by default.
+    Made building on homebrew installations a bit better.
+    Remove temporary variables for glog and gflags.
+
+commit eaccbb345614c0d24c5e21fa931f470cfda874df
+Author: Keir Mierle <mierle@gmail.com>
+Date:   Wed May 9 05:31:29 2012 -0700
+
+    Remove unused template parameter from VariadicEvaluate.
+
+commit 82f4b88c34b0b2cf85064e5fc20e374e978b2e3b
+Author: Sameer Agarwal <sameeragarwal@google.com>
+Date:   Sun May 6 21:05:28 2012 -0700
+
+    Extend support writing linear least squares problems to disk.
+    
+    1. Make the mechanism for writing problems to disk, generic and
+    controllable using an enum DumpType visible in the API.
+    
+    2. Instead of single file containing protocol buffers, now matrices can
+    be written in a matlab/octave friendly format. This is now the default.
+    
+    3. The support for writing problems to disk is moved into
+    linear_least_squares_problem.cc/h
+    
+    4. SparseMatrix now has a ToTextFile virtual method which is
+    implemented by each of its subclasses to write a (i,j,s) triplets.
+    
+    5. Minor changes to simple_bundle_adjuster to enable logging at startup.
+
 commit d297f8d3d3f5025c24752f0f4c1ec2469a769f99
 Merge: 7e74d81 f8bd7fa
 Author: Keir Mierle <mierle@gmail.com>
index 0ac6240dfabd278e42d9d7c45e8aa893526f5a35..da9ee2c7993e7b5b4816838d0aca32368ec3fdd1 100644 (file)
 // "MyScalarCostFunction", "1, 2, 2", describe the functor as computing a
 // 1-dimensional output from two arguments, both 2-dimensional.
 //
+// The autodiff cost function also supports cost functions with a
+// runtime-determined number of residuals. For example:
+//
+//   CostFunction* cost_function
+//       = new AutoDiffCostFunction<MyScalarCostFunction, DYNAMIC, 2, 2>(
+//           new CostFunctionWithDynamicNumResiduals(1.0),   ^     ^  ^
+//           runtime_number_of_residuals); <----+            |     |  |
+//                                              |            |     |  |
+//                                              |            |     |  |
+//             Actual number of residuals ------+            |     |  |
+//             Indicate dynamic number of residuals ---------+     |  |
+//             Dimension of x -------------------------------------+  |
+//             Dimension of y ----------------------------------------+
+//
 // The framework can currently accommodate cost functions of up to 6 independent
 // variables, and there is no limit on the dimensionality of each of them.
 //
 #include "ceres/internal/autodiff.h"
 #include "ceres/internal/scoped_ptr.h"
 #include "ceres/sized_cost_function.h"
+#include "ceres/types.h"
 
 namespace ceres {
 
-// A cost function which computes the derivative of the cost with respect to the
-// parameters (a.k.a. the jacobian) using an autodifferentiation framework. The
-// first template argument is the functor object, described in the header
-// comment. The second argument is the dimension of the residual, and subsequent
+// A cost function which computes the derivative of the cost with respect to
+// the parameters (a.k.a. the jacobian) using an autodifferentiation framework.
+// The first template argument is the functor object, described in the header
+// comment. The second argument is the dimension of the residual (or
+// ceres::DYNAMIC to indicate it will be set at runtime), and subsequent
 // arguments describe the size of the Nth parameter, one per parameter.
 //
-// The constructor, which takes a cost functor, takes ownership of the functor.
+// The constructors take ownership of the cost functor.
+//
+// If the number of residuals (argument "M" below) is ceres::DYNAMIC, then the
+// two-argument constructor must be used. The second constructor takes a number
+// of residuals (in addition to the templated number of residuals). This allows
+// for varying the number of residuals for a single autodiff cost function at
+// runtime.
 template <typename CostFunctor,
-          int M,        // Number of residuals.
+          int M,        // Number of residuals, or ceres::DYNAMIC.
           int N0,       // Number of parameters in block 0.
           int N1 = 0,   // Number of parameters in block 1.
           int N2 = 0,   // Number of parameters in block 2.
@@ -136,8 +158,25 @@ template <typename CostFunctor,
 class AutoDiffCostFunction :
   public SizedCostFunction<M, N0, N1, N2, N3, N4, N5> {
  public:
-  // Takes ownership of functor.
-  explicit AutoDiffCostFunction(CostFunctor* functor) : functor_(functor) {}
+  // Takes ownership of functor. Uses the template-provided value for the
+  // number of residuals ("M").
+  explicit AutoDiffCostFunction(CostFunctor* functor)
+      : functor_(functor) {
+    CHECK_NE(M, DYNAMIC) << "Can't run the fixed-size constructor if the "
+                          << "number of residuals is set to ceres::DYNAMIC.";
+  }
+
+  // Takes ownership of functor. Ignores the template-provided number of
+  // residuals ("M") in favor of the "num_residuals" argument provided.
+  //
+  // This allows for having autodiff cost functions which return varying
+  // numbers of residuals at runtime.
+  AutoDiffCostFunction(CostFunctor* functor, int num_residuals)
+      : functor_(functor) {
+    CHECK_EQ(M, DYNAMIC) << "Can't run the dynamic-size constructor if the "
+                          << "number of residuals is not ceres::DYNAMIC.";
+    SizedCostFunction<M, N0, N1, N2, N3, N4, N5>::set_num_residuals(num_residuals);
+  }
 
   virtual ~AutoDiffCostFunction() {}
 
@@ -151,14 +190,16 @@ class AutoDiffCostFunction :
                         double** jacobians) const {
     if (!jacobians) {
       return internal::VariadicEvaluate<
-          CostFunctor, double, M, N0, N1, N2, N3, N4, N5>
+          CostFunctor, double, N0, N1, N2, N3, N4, N5>
           ::Call(*functor_, parameters, residuals);
     }
     return internal::AutoDiff<CostFunctor, double,
-           M, N0, N1, N2, N3, N4, N5>::Differentiate(*functor_,
-                                                     parameters,
-                                                     residuals,
-                                                     jacobians);
+           N0, N1, N2, N3, N4, N5>::Differentiate(
+               *functor_,
+               parameters,
+               SizedCostFunction<M, N0, N1, N2, N3, N4, N5>::num_residuals(),
+               residuals,
+               jacobians);
   }
 
  private:
index 1a9d396c9efc837f2f09a2e339dd3825468c9f8a..4f5081f8f66327318289155a8179d1ac657291d5 100644 (file)
 
 #include <glog/logging.h>
 #include "ceres/jet.h"
+#include "ceres/internal/eigen.h"
 #include "ceres/internal/fixed_array.h"
 
 namespace ceres {
@@ -185,18 +186,12 @@ inline void Take0thOrderPart(int M, const JetT *src, T dst) {
 
 // Takes N 1st order parts, starting at index N0, and puts them in the M x N
 // matrix 'dst'. This is used to pick out the "matrix" parts of the extended y.
-template <typename JetT, typename T, int M, int N0, int N>
-inline void Take1stOrderPart(const JetT *src, T *dst) {
+template <typename JetT, typename T, int N0, int N>
+inline void Take1stOrderPart(const int M, const JetT *src, T *dst) {
   DCHECK(src);
   DCHECK(dst);
-  // TODO(keir): Change Jet to use a single array, where v[0] is the
-  // non-infinitesimal part rather than "a". That way it's possible to use a
-  // single memcpy or eigen operation, rather than the explicit loop. The loop
-  // doesn't exploit any SSE or other intrinsics.
   for (int i = 0; i < M; ++i) {
-    for (int j = 0; j < N; ++j) {
-      dst[N * i + j] = src[i].v[N0 + j];
-    }
+    Eigen::Map<Eigen::Matrix<T, N, 1> >(dst + N * i, N) = src[i].v.template segment<N>(N0);
   }
 }
 
@@ -208,7 +203,7 @@ inline void Take1stOrderPart(const JetT *src, T *dst) {
 // Supporting variadic functions is the primary source of complexity in the
 // autodiff implementation.
 
-template<typename Functor, typename T, int kNumOutputs,
+template<typename Functor, typename T,
          int N0, int N1, int N2, int N3, int N4, int N5>
 struct VariadicEvaluate {
   static bool Call(const Functor& functor, T const *const *input, T* output) {
@@ -222,9 +217,9 @@ struct VariadicEvaluate {
   }
 };
 
-template<typename Functor, typename T, int kNumOutputs,
+template<typename Functor, typename T,
          int N0, int N1, int N2, int N3, int N4>
-struct VariadicEvaluate<Functor, T, kNumOutputs, N0, N1, N2, N3, N4, 0> {
+struct VariadicEvaluate<Functor, T, N0, N1, N2, N3, N4, 0> {
   static bool Call(const Functor& functor, T const *const *input, T* output) {
     return functor(input[0],
                    input[1],
@@ -235,9 +230,9 @@ struct VariadicEvaluate<Functor, T, kNumOutputs, N0, N1, N2, N3, N4, 0> {
   }
 };
 
-template<typename Functor, typename T, int kNumOutputs,
+template<typename Functor, typename T,
          int N0, int N1, int N2, int N3>
-struct VariadicEvaluate<Functor, T, kNumOutputs, N0, N1, N2, N3, 0, 0> {
+struct VariadicEvaluate<Functor, T, N0, N1, N2, N3, 0, 0> {
   static bool Call(const Functor& functor, T const *const *input, T* output) {
     return functor(input[0],
                    input[1],
@@ -247,9 +242,9 @@ struct VariadicEvaluate<Functor, T, kNumOutputs, N0, N1, N2, N3, 0, 0> {
   }
 };
 
-template<typename Functor, typename T, int kNumOutputs,
+template<typename Functor, typename T,
          int N0, int N1, int N2>
-struct VariadicEvaluate<Functor, T, kNumOutputs, N0, N1, N2, 0, 0, 0> {
+struct VariadicEvaluate<Functor, T, N0, N1, N2, 0, 0, 0> {
   static bool Call(const Functor& functor, T const *const *input, T* output) {
     return functor(input[0],
                    input[1],
@@ -258,9 +253,9 @@ struct VariadicEvaluate<Functor, T, kNumOutputs, N0, N1, N2, 0, 0, 0> {
   }
 };
 
-template<typename Functor, typename T, int kNumOutputs,
+template<typename Functor, typename T,
          int N0, int N1>
-struct VariadicEvaluate<Functor, T, kNumOutputs, N0, N1, 0, 0, 0, 0> {
+struct VariadicEvaluate<Functor, T, N0, N1, 0, 0, 0, 0> {
   static bool Call(const Functor& functor, T const *const *input, T* output) {
     return functor(input[0],
                    input[1],
@@ -268,8 +263,8 @@ struct VariadicEvaluate<Functor, T, kNumOutputs, N0, N1, 0, 0, 0, 0> {
   }
 };
 
-template<typename Functor, typename T, int kNumOutputs, int N0>
-struct VariadicEvaluate<Functor, T, kNumOutputs, N0, 0, 0, 0, 0, 0> {
+template<typename Functor, typename T, int N0>
+struct VariadicEvaluate<Functor, T, N0, 0, 0, 0, 0, 0> {
   static bool Call(const Functor& functor, T const *const *input, T* output) {
     return functor(input[0],
                    output);
@@ -279,11 +274,12 @@ struct VariadicEvaluate<Functor, T, kNumOutputs, N0, 0, 0, 0, 0, 0> {
 // This is in a struct because default template parameters on a function are not
 // supported in C++03 (though it is available in C++0x). N0 through N5 are the
 // dimension of the input arguments to the user supplied functor.
-template <typename Functor, typename T, int kNumOutputs,
+template <typename Functor, typename T,
           int N0 = 0, int N1 = 0, int N2 = 0, int N3 = 0, int N4 = 0, int N5=0>
 struct AutoDiff {
   static bool Differentiate(const Functor& functor,
                             T const *const *parameters,
+                            int num_outputs,
                             T *function_value,
                             T **jacobians) {
     typedef Jet<T, N0 + N1 + N2 + N3 + N4 + N5> JetT;
@@ -300,10 +296,10 @@ struct AutoDiff {
         << "(ignore trailing 0s): " << N0 << ", " << N1 << ", " << N2 << ", "
         << N3 << ", " << N4 << ", " << N5;
 
-    DCHECK_GT(kNumOutputs, 0);
+    DCHECK_GT(num_outputs, 0);
 
     FixedArray<JetT, (256 * 7) / sizeof(JetT)> x(
-        N0 + N1 + N2 + N3 + N4 + N5 + kNumOutputs);
+        N0 + N1 + N2 + N3 + N4 + N5 + num_outputs);
 
     // It's ugly, but it works.
     const int jet0 = 0;
@@ -339,22 +335,22 @@ struct AutoDiff {
     CERES_MAKE_1ST_ORDER_PERTURBATION(5);
 #undef CERES_MAKE_1ST_ORDER_PERTURBATION
 
-    if (!VariadicEvaluate<Functor, JetT, kNumOutputs,
+    if (!VariadicEvaluate<Functor, JetT,
                           N0, N1, N2, N3, N4, N5>::Call(
         functor, unpacked_parameters, output)) {
       return false;
     }
 
-    internal::Take0thOrderPart(kNumOutputs, output, function_value);
+    internal::Take0thOrderPart(num_outputs, output, function_value);
 
 #define CERES_TAKE_1ST_ORDER_PERTURBATION(i) \
     if (N ## i) { \
       if (jacobians[i]) { \
         internal::Take1stOrderPart<JetT, T, \
-                                   kNumOutputs, \
                                    jet ## i, \
-                                   N ## i>(output, \
-                                          jacobians[i]); \
+                                   N ## i>(num_outputs, \
+                                           output, \
+                                           jacobians[i]); \
       } \
     }
     CERES_TAKE_1ST_ORDER_PERTURBATION(0);
index 968285b8f1e05ad45b2b8eff025425f2e6ea5913..2894a9fba5c0607a279eb89416640635333e1f1f 100644 (file)
 //
 // A convenience class for cost functions which are statically sized.
 // Compared to the dynamically-sized base class, this reduces boilerplate.
+//
+// The kNumResiduals template parameter can be a constant such as 2 or 5, or it
+// can be ceres::DYNAMIC. If kNumResiduals is ceres::DYNAMIC, then subclasses
+// are responsible for calling set_num_residuals() at runtime.
 
 #ifndef CERES_PUBLIC_SIZED_COST_FUNCTION_H_
 #define CERES_PUBLIC_SIZED_COST_FUNCTION_H_
 
 #include <glog/logging.h>
+#include "ceres/types.h"
 #include "ceres/cost_function.h"
 
 namespace ceres {
@@ -45,11 +50,12 @@ class SizedCostFunction : public CostFunction {
  public:
   SizedCostFunction() {
     // Sanity checking.
-    DCHECK_GT(kNumResiduals, 0) << "Cost functions must have at least "
-                                << "one residual block.";
-    DCHECK_GT(N0, 0)
+    CHECK(kNumResiduals > 0 || kNumResiduals == DYNAMIC)
+        << "Cost functions must have at least one residual block.";
+
+    CHECK_GT(N0, 0)
         << "Cost functions must have at least one parameter block.";
-    DCHECK((!N1 && !N2 && !N3 && !N4 && !N5) ||
+    CHECK((!N1 && !N2 && !N3 && !N4 && !N5) ||
            ((N1 > 0) && !N2 && !N3 && !N4 && !N5) ||
            ((N1 > 0) && (N2 > 0) && !N3 && !N4 && !N5) ||
            ((N1 > 0) && (N2 > 0) && (N3 > 0) && !N4 && !N5) ||
index 15fd7332d2122dde9b95f0e4eaa0f005b8aad8a4..bd6692720234bf670ff2f9efae67fe5bddde02c7 100644 (file)
@@ -83,7 +83,8 @@ class Solver {
       minimizer_progress_to_stdout = false;
       return_initial_residuals = false;
       return_final_residuals = false;
-      lsqp_dump_format = "lm_iteration_%03d.lsqp";
+      lsqp_dump_directory = "/tmp";
+      lsqp_dump_format_type = TEXTFILE;
       crash_and_dump_lsqp_on_failure = false;
       check_gradients = false;
       gradient_check_relative_precision = 1e-8;
@@ -213,12 +214,8 @@ class Solver {
     //
     // This is ignored if protocol buffers are disabled.
     vector<int> lsqp_iterations_to_dump;
-
-    // Format string for the file name used for dumping the least
-    // squares problem to disk. If the format is 'ascii', then the
-    // problem is logged to the screen; don't try this with large
-    // problems or expect a frozen terminal.
-    string lsqp_dump_format;
+    string lsqp_dump_directory;
+    DumpFormatType lsqp_dump_format_type;
 
     // Dump the linear least squares problem to disk if the minimizer
     // fails due to NUMERICAL_FAILURE and crash the process. This flag
index b83a266d1ba9e6c66544fae5d098ff6dc9d10789..a30c79029ace88033244c0516632db9aec183db7 100644 (file)
@@ -210,6 +210,40 @@ enum CallbackReturnType {
   SOLVER_TERMINATE_SUCCESSFULLY
 };
 
+// The format in which linear least squares problems should be logged
+// when Solver::Options::lsqp_iterations_to_dump is non-empty.
+enum DumpFormatType {
+  // Print the linear least squares problem in a human readable format
+  // to stderr. The Jacobian is printed as a dense matrix. The vectors
+  // D, x and f are printed as dense vectors. This should only be used
+  // for small problems.
+  CONSOLE,
+
+  // Write out the linear least squares problem to the directory
+  // pointed to by Solver::Options::lsqp_dump_directory as a protocol
+  // buffer. linear_least_squares_problems.h/cc contains routines for
+  // loading these problems. For details on the on disk format used,
+  // see matrix.proto. The files are named lm_iteration_???.lsqp.
+  PROTOBUF,
+
+  // Write out the linear least squares problem to the directory
+  // pointed to by Solver::Options::lsqp_dump_directory as text files
+  // which can be read into MATLAB/Octave. The Jacobian is dumped as a
+  // text file containing (i,j,s) triplets, the vectors D, x and f are
+  // dumped as text files containing a list of their values.
+  //
+  // A MATLAB/octave script called lm_iteration_???.m is also output,
+  // which can be used to parse and load the problem into memory.
+  TEXTFILE
+};
+
+// For SizedCostFunction and AutoDiffCostFunction, DYNAMIC can be specified for
+// the number of residuals. If specified, then the number of residuas for that
+// cost function can vary at runtime.
+enum DimensionType {
+  DYNAMIC = -1
+};
+
 const char* LinearSolverTypeToString(LinearSolverType type);
 const char* PreconditionerTypeToString(PreconditionerType type);
 const char* LinearSolverTerminationTypeToString(
index c1be9402b7863a44feda7fbd267233321b46237d..7dd395e29751b952c47fe6072cf6d68765044eb6 100644 (file)
@@ -259,5 +259,28 @@ void BlockSparseMatrix::ToProto(SparseMatrixProto* outer_proto) const {
 }
 #endif
 
+void BlockSparseMatrix::ToTextFile(FILE* file) const {
+  CHECK_NOTNULL(file);
+  for (int i = 0; i < block_structure_->rows.size(); ++i) {
+    const int row_block_pos = block_structure_->rows[i].block.position;
+    const int row_block_size = block_structure_->rows[i].block.size;
+    const vector<Cell>& cells = block_structure_->rows[i].cells;
+    for (int j = 0; j < cells.size(); ++j) {
+      const int col_block_id = cells[j].block_id;
+      const int col_block_size = block_structure_->cols[col_block_id].size;
+      const int col_block_pos = block_structure_->cols[col_block_id].position;
+      int jac_pos = cells[j].position;
+      for (int r = 0; r < row_block_size; ++r) {
+        for (int c = 0; c < col_block_size; ++c) {
+          fprintf(file, "% 10d % 10d %17f\n",
+                  row_block_pos + r,
+                  col_block_pos + c,
+                  values_[jac_pos++]);
+        }
+      }
+    }
+  }
+}
+
 }  // namespace internal
 }  // namespace ceres
index b151dd0e248a0ef5f606bba429c65f003c726569..f71446e8f58ac3d24095c3ecdd44479b821ca73c 100644 (file)
@@ -113,6 +113,8 @@ class BlockSparseMatrix : public BlockSparseMatrixBase {
 #ifndef CERES_DONT_HAVE_PROTOCOL_BUFFERS
   virtual void ToProto(SparseMatrixProto* proto) const;
 #endif
+  virtual void ToTextFile(FILE* file) const;
+
   virtual int num_rows()         const { return num_rows_;     }
   virtual int num_cols()         const { return num_cols_;     }
   virtual int num_nonzeros()     const { return num_nonzeros_; }
index 8fd568ffcc3db6977afa2f660992a00fc78e5a0b..95edf5396af83fdf91b95b1c64a78e61c661202c 100644 (file)
@@ -321,5 +321,14 @@ void CompressedRowSparseMatrix::AppendRows(const CompressedRowSparseMatrix& m) {
   num_rows_ += m.num_rows();
 }
 
+void CompressedRowSparseMatrix::ToTextFile(FILE* file) const {
+  CHECK_NOTNULL(file);
+  for (int r = 0; r < num_rows_; ++r) {
+    for (int idx = rows_[r]; idx < rows_[r + 1]; ++idx) {
+      fprintf(file, "% 10d % 10d %17f\n", r, cols_[idx], values_[idx]);
+    }
+  }
+}
+
 }  // namespace internal
 }  // namespace ceres
index 43712a866405715549d389398a51a8dec7492c52..9a39d28e111b0d69cde03c6a231eac0fc6393052 100644 (file)
@@ -88,7 +88,7 @@ class CompressedRowSparseMatrix : public SparseMatrix {
 #ifndef CERES_DONT_HAVE_PROTOCOL_BUFFERS
   virtual void ToProto(SparseMatrixProto* proto) const;
 #endif
-
+  virtual void ToTextFile(FILE* file) const;
   virtual int num_rows() const { return num_rows_; }
   virtual int num_cols() const { return num_cols_; }
   virtual int num_nonzeros() const { return rows_[num_rows_]; }
index ffbfab61de1917611286d6284ceb59f3b1fe1c36..5d392ba6c3bf14fc2a407ffe183b6f2173e09ee4 100644 (file)
@@ -179,6 +179,19 @@ AlignedMatrixRef DenseSparseMatrix::mutable_matrix() {
   return AlignedMatrixRef(m_.data(), m_.rows(), m_.cols());
 }
 
+void DenseSparseMatrix::ToTextFile(FILE* file) const {
+  CHECK_NOTNULL(file);
+  const int active_rows =
+      (has_diagonal_reserved_ && !has_diagonal_appended_)
+      ? (m_.rows() - m_.cols()) 
+      : m_.rows();
+
+  for (int r = 0; r < active_rows; ++r) {
+    for (int c = 0; c < m_.cols(); ++c) {
+      fprintf(file,  "% 10d % 10d %17f\n", r, c, m_(r, c));
+    }
+  }
+}
 
 }  // namespace internal
 }  // namespace ceres
index 5ce29eef51bc32e1221604af3ba4aa2bc89b7260..416c2143c2c774a7042b2efc9442ceef935a4d2c 100644 (file)
@@ -70,6 +70,7 @@ class DenseSparseMatrix : public SparseMatrix {
 #ifndef CERES_DONT_HAVE_PROTOCOL_BUFFERS
   virtual void ToProto(SparseMatrixProto* proto) const;
 #endif
+  virtual void ToTextFile(FILE* file) const;
   virtual int num_rows() const;
   virtual int num_cols() const;
   virtual int num_nonzeros() const;
index 3ad359e63e4fc327606ab73c18ce0093119d541c..b40a5162adcae0c4ffb3321c5cef5a36b760c73d 100644 (file)
@@ -58,6 +58,7 @@
 #include "Eigen/Core"
 #include "ceres/evaluator.h"
 #include "ceres/file.h"
+#include "ceres/linear_least_squares_problems.h"
 #include "ceres/linear_solver.h"
 #include "ceres/matrix_proto.h"
 #include "ceres/sparse_matrix.h"
@@ -107,62 +108,6 @@ void LevenbergMarquardtDiagonal(const SparseMatrix& jacobian,
   }
 }
 
-string DumpLinearSolverProblem(
-    int iteration,
-    const SparseMatrix* A,
-    const double* D,
-    const double* b,
-    const double* x,
-    const Minimizer::Options& solver_options) {
-  if (solver_options.lsqp_dump_format == "ascii") {
-    // Dump to the screen instead of to file. Useful for debugging.
-    Matrix AA;
-    A->ToDenseMatrix(&AA);
-    LOG(INFO) << "A^T: \n" << AA.transpose();
-    if (D) {
-      LOG(INFO) << "A's appended diagonal:\n"
-                << ConstVectorRef(D, A->num_cols());
-    }
-    LOG(INFO) << "b: \n" << ConstVectorRef(b, A->num_rows());
-    LOG(INFO) << "x: \n" << ConstVectorRef(x, A->num_cols());
-    return "";
-  }
-#ifndef CERES_DONT_HAVE_PROTOCOL_BUFFERS
-  LinearLeastSquaresProblemProto lsqp;
-  A->ToProto(lsqp.mutable_a());
-  for (int i = 0; i < A->num_rows(); ++i) {
-    lsqp.add_b(b[i]);
-  }
-  if (D) {
-    for (int i = 0; i < A->num_cols(); ++i) {
-      lsqp.add_d(D[i]);
-    }
-  }
-  if (x) {
-    for (int i = 0; i < A->num_cols(); ++i) {
-      lsqp.add_x(x[i]);
-    }
-  }
-
-  lsqp.set_num_eliminate_blocks(solver_options.num_eliminate_blocks);
-
-  CHECK(solver_options.lsqp_dump_format.size());
-  string filename =
-      StringPrintf(solver_options.lsqp_dump_format.c_str(),  // NOLINT
-                   iteration);
-  VLOG(1) << "Dumping least squares problem for iteration " << iteration
-          << " to disk. File: " << filename;
-  WriteStringToFileOrDie(lsqp.SerializeAsString(), filename);
-  VLOG(2) << "Done dumping to disk";
-  return filename;
-#else
-  LOG(ERROR) << "Dumping least squares problems is only "
-             << "supported when Ceres is compiled with "
-             << "protocol buffer support.";
-  return "";
-#endif
-}
-
 bool RunCallback(IterationCallback* callback,
                  const IterationSummary& iteration_summary,
                  Solver::Summary* summary) {
@@ -380,12 +325,17 @@ void LevenbergMarquardt::Minimize(const Minimizer::Options& options,
       if (binary_search(iterations_to_dump.begin(),
                         iterations_to_dump.end(),
                         iteration)) {
-        DumpLinearSolverProblem(iteration,
-                                jacobian.get(),
-                                muD.data(),
-                                f.data(),
-                                lm_step.data(),
-                                options);
+        CHECK(DumpLinearLeastSquaresProblem(options.lsqp_dump_directory,
+                                            iteration,
+                                            options.lsqp_dump_format_type,
+                                            jacobian.get(),
+                                            muD.data(),
+                                            f.data(),
+                                            lm_step.data(),
+                                            options.num_eliminate_blocks))
+            << "Tried writing linear least squares problem: " 
+            << options.lsqp_dump_directory
+            << " but failed.";
       }
 
       // We ignore the case where the linear solver did not converge,
@@ -413,7 +363,7 @@ void LevenbergMarquardt::Minimize(const Minimizer::Options& options,
           (x_norm + options.parameter_tolerance)) {
         summary->termination_type = PARAMETER_TOLERANCE;
         VLOG(1) << "Terminating on PARAMETER_TOLERANCE."
-             << "Relative step size: " << step_norm / step_size_tolerance
+                << "Relative step size: " << step_norm / step_size_tolerance
             << " <= " << options.parameter_tolerance;
         return;
       }
@@ -545,33 +495,37 @@ void LevenbergMarquardt::Minimize(const Minimizer::Options& options,
     }
 
     if (num_consecutive_insane_steps == kMaxLinearSolverRetries) {
-      VLOG(1) << "Too many consecutive retries; ending with numerical fail.";
       summary->termination_type = NUMERICAL_FAILURE;
+      VLOG(1) << "Too many consecutive retries; ending with numerical fail.";
 
       if (!options.crash_and_dump_lsqp_on_failure) {
         return;
       }
 
       // Dump debugging information to disk.
-      CHECK(!options.lsqp_dump_format.empty())
+      CHECK(options.lsqp_dump_format_type == TEXTFILE ||
+            options.lsqp_dump_format_type == PROTOBUF)
           << "Dumping the linear least squares problem on crash "
-          << "requires Solver::Options::lsqp_dump_format set a "
-          << "filename";
-      CHECK_NE(options.lsqp_dump_format, "ascii")
-          << "Dumping the linear least squares problem on crash "
-          << "requires Solver::Options::lsqp_dump_format set a "
-          << "filename";
-
-      const string filename = DumpLinearSolverProblem(iteration,
-                                                      jacobian.get(),
-                                                      muD.data(),
-                                                      f.data(),
-                                                      lm_step.data(),
-                                                      options);
-      LOG(FATAL) << "Linear least squares problem saved to " << filename
-                 << " please provide this to the Ceres developers for "
-                 << " debugging along with the v=2 log.";
-      return;
+          << "requires Solver::Options::lsqp_dump_format_type to be "
+          << "PROTOBUF or TEXTFILE.";
+
+      if (DumpLinearLeastSquaresProblem(options.lsqp_dump_directory,
+                                        iteration,
+                                        options.lsqp_dump_format_type,
+                                        jacobian.get(),
+                                        muD.data(),
+                                        f.data(),
+                                        lm_step.data(),
+                                        options.num_eliminate_blocks)) {
+        LOG(FATAL) << "Linear least squares problem saved to: " 
+                   << options.lsqp_dump_directory
+                   << ". Please provide this to the Ceres developers for "
+                   << " debugging along with the v=2 log.";
+      } else {
+        LOG(FATAL) << "Tried writing linear least squares problem: " 
+                   << options.lsqp_dump_directory
+                   << " but failed.";
+      }
     }
 
     if (!step_is_successful) {
index 9fc5ff8a1c7dcbec1590b5f9dd541fd8f8038664..cca9f442fe711ef71a0a872c2b13649941c3a606 100644 (file)
 
 #include "ceres/linear_least_squares_problems.h"
 
+#include <cstdio>
 #include <string>
 #include <vector>
 #include <glog/logging.h>
 #include "ceres/block_sparse_matrix.h"
 #include "ceres/block_structure.h"
+#include "ceres/casts.h"
 #include "ceres/compressed_row_sparse_matrix.h"
 #include "ceres/file.h"
 #include "ceres/matrix_proto.h"
 #include "ceres/triplet_sparse_matrix.h"
+#include "ceres/stringprintf.h"
 #include "ceres/internal/scoped_ptr.h"
 #include "ceres/types.h"
 
@@ -570,5 +573,172 @@ LinearLeastSquaresProblem* LinearLeastSquaresProblem3() {
   return problem;
 }
 
+bool DumpLinearLeastSquaresProblemToConsole(const string& directory,
+                                            int iteration,
+                                            const SparseMatrix* A,
+                                            const double* D,
+                                            const double* b,
+                                            const double* x,
+                                            int num_eliminate_blocks) {
+  CHECK_NOTNULL(A);
+  Matrix AA;
+  A->ToDenseMatrix(&AA);
+  LOG(INFO) << "A^T: \n" << AA.transpose();
+
+  if (D != NULL) {
+    LOG(INFO) << "A's appended diagonal:\n"
+              << ConstVectorRef(D, A->num_cols());
+  }
+
+  if (b != NULL) {
+    LOG(INFO) << "b: \n" << ConstVectorRef(b, A->num_rows());
+  }
+
+  if (x != NULL) {
+    LOG(INFO) << "x: \n" << ConstVectorRef(x, A->num_cols());
+  }
+  return true;
+};
+
+#ifndef CERES_DONT_HAVE_PROTOCOL_BUFFERS
+bool DumpLinearLeastSquaresProblemToProtocolBuffer(const string& directory,
+                                                   int iteration,
+                                                   const SparseMatrix* A,
+                                                   const double* D,
+                                                   const double* b,
+                                                   const double* x,
+                                                   int num_eliminate_blocks) {
+  CHECK_NOTNULL(A);
+  LinearLeastSquaresProblemProto lsqp;
+  A->ToProto(lsqp.mutable_a());
+
+  if (D != NULL) {
+    for (int i = 0; i < A->num_cols(); ++i) {
+      lsqp.add_d(D[i]);
+    }
+  }
+
+  if (b != NULL) {
+    for (int i = 0; i < A->num_rows(); ++i) {
+      lsqp.add_b(b[i]);
+    }
+  }
+
+  if (x != NULL) {
+    for (int i = 0; i < A->num_cols(); ++i) {
+      lsqp.add_x(x[i]);
+    }
+  }
+
+  lsqp.set_num_eliminate_blocks(num_eliminate_blocks);
+  string format_string = JoinPath(directory,
+                                  "lm_iteration_%03d.lsqp");
+  string filename =
+      StringPrintf(format_string.c_str(),  iteration);
+  LOG(INFO) << "Dumping least squares problem for iteration " << iteration
+            << " to disk. File: " << filename;
+  WriteStringToFileOrDie(lsqp.SerializeAsString(), filename);
+  return true;
+}
+#else
+bool DumpLinearLeastSquaresProblemToProtocolBuffer(const string& directory,
+                                                   int iteration,
+                                                   const SparseMatrix* A,
+                                                   const double* D,
+                                                   const double* b,
+                                                   const double* x,
+                                                   int num_eliminate_blocks) {
+  LOG(ERROR) << "Dumping least squares problems is only "
+             << "supported when Ceres is compiled with "
+             << "protocol buffer support.";
+  return false;
+}
+#endif
+
+void WriteArrayToFileOrDie(const string& filename,
+                           const double* x,
+                           const int size) {
+  CHECK_NOTNULL(x);
+  VLOG(2) << "Writing array to: " << filename;
+  FILE* fptr = fopen(filename.c_str(), "w");
+  CHECK_NOTNULL(fptr);
+  for (int i = 0; i < size; ++i) {
+    fprintf(fptr, "%17f\n", x[i]);
+  }
+  fclose(fptr);
+}
+
+bool DumpLinearLeastSquaresProblemToTextFile(const string& directory,
+                                             int iteration,
+                                             const SparseMatrix* A,
+                                             const double* D,
+                                             const double* b,
+                                             const double* x,
+                                             int num_eliminate_blocks) {
+  CHECK_NOTNULL(A);
+  string format_string = JoinPath(directory,
+                                  "lm_iteration_%03d");
+  string filename_prefix =
+      StringPrintf(format_string.c_str(), iteration);
+
+  {
+    string filename = filename_prefix + "_A.txt";
+    LOG(INFO) << "writing to: " << filename;
+    FILE* fptr = fopen(filename.c_str(), "w");
+    CHECK_NOTNULL(fptr);
+    A->ToTextFile(fptr);
+    fclose(fptr);
+  }
+
+  if (D != NULL) {
+    string filename = filename_prefix + "_D.txt";
+    WriteArrayToFileOrDie(filename, D, A->num_cols());
+  }
+
+  if (b != NULL) {
+    string filename = filename_prefix + "_b.txt";
+    WriteArrayToFileOrDie(filename, b, A->num_rows());
+  }
+
+  if (x != NULL) {
+    string filename = filename_prefix + "_x.txt";
+    WriteArrayToFileOrDie(filename, x, A->num_cols());
+  }
+
+  return true;
+}
+
+bool DumpLinearLeastSquaresProblem(const string& directory,
+                                  int iteration,
+                                   DumpFormatType dump_format_type,
+                                   const SparseMatrix* A,
+                                   const double* D,
+                                   const double* b,
+                                   const double* x,
+                                   int num_eliminate_blocks) {
+  switch (dump_format_type) {
+    case (CONSOLE):
+      return DumpLinearLeastSquaresProblemToConsole(directory,
+                                                    iteration,
+                                                    A, D, b, x,
+                                                    num_eliminate_blocks);
+    case (PROTOBUF):
+      return DumpLinearLeastSquaresProblemToProtocolBuffer(
+          directory,
+          iteration,
+          A, D, b, x,
+          num_eliminate_blocks);
+    case (TEXTFILE):
+      return DumpLinearLeastSquaresProblemToTextFile(directory,
+                                                     iteration,
+                                                     A, D, b, x,
+                                                     num_eliminate_blocks);
+    default:
+      LOG(FATAL) << "Unknown DumpFormatType " << dump_format_type;
+  };
+
+  return true;
+}
+
 }  // namespace internal
 }  // namespace ceres
index 46a624bd73fc15e7133125884fb0ab2bd32ea7fa..553cc0d3db34dd3f8ac3ad76c8c32c9399a8de58 100644 (file)
@@ -32,7 +32,7 @@
 #define CERES_INTERNAL_LINEAR_LEAST_SQUARES_PROBLEMS_H_
 
 #include <string>
-
+#include <vector>
 #include "ceres/sparse_matrix.h"
 #include "ceres/internal/port.h"
 #include "ceres/internal/scoped_ptr.h"
@@ -71,6 +71,16 @@ LinearLeastSquaresProblem* LinearLeastSquaresProblem1();
 LinearLeastSquaresProblem* LinearLeastSquaresProblem2();
 LinearLeastSquaresProblem* LinearLeastSquaresProblem3();
 
+// Write the linear least squares problem to disk. The exact format
+// depends on dump_format_type.
+bool DumpLinearLeastSquaresProblem(const string& directory,
+                                  int iteration,
+                                   DumpFormatType dump_format_type,
+                                   const SparseMatrix* A,
+                                   const double* D,
+                                   const double* b,
+                                   const double* x,
+                                   int num_eliminate_blocks);
 }  // namespace internal
 }  // namespace ceres
 
index 71163a8ea6f6bbc6d331cb7787d6537ef6c940ad..77cb00cb6b4874f8352e8898fdae5780fac69520 100644 (file)
@@ -59,8 +59,9 @@ class Minimizer {
       tau = options.tau;
       jacobi_scaling = options.jacobi_scaling;
       crash_and_dump_lsqp_on_failure = options.crash_and_dump_lsqp_on_failure;
-      lsqp_dump_format = options.lsqp_dump_format;
+      lsqp_dump_directory = options.lsqp_dump_directory;
       lsqp_iterations_to_dump = options.lsqp_iterations_to_dump;
+      lsqp_dump_format_type = options.lsqp_dump_format_type;
       num_eliminate_blocks = options.num_eliminate_blocks;
       logging_type = options.logging_type;
     }
@@ -75,8 +76,9 @@ class Minimizer {
     double tau;
     bool jacobi_scaling;
     bool crash_and_dump_lsqp_on_failure;
-    string lsqp_dump_format;
     vector<int> lsqp_iterations_to_dump;
+    DumpFormatType lsqp_dump_format_type;
+    string lsqp_dump_directory;
     int num_eliminate_blocks;
     LoggingType logging_type;
 
index 962b803dd877cd0bddffdd1c70513d9fa14e732b..562210dfec8ccb4a651de62bd89edcadb1589cf6 100644 (file)
@@ -33,6 +33,7 @@
 #ifndef CERES_INTERNAL_SPARSE_MATRIX_H_
 #define CERES_INTERNAL_SPARSE_MATRIX_H_
 
+#include <cstdio>
 #include "ceres/linear_operator.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/types.h"
@@ -87,9 +88,14 @@ class SparseMatrix : public LinearOperator {
 
 #ifndef CERES_DONT_HAVE_PROTOCOL_BUFFERS
   // Dump the sparse matrix to a proto. Destroys the contents of proto.
-  virtual void ToProto(SparseMatrixProto *proto) const = 0;
+  virtual void ToProto(SparseMatrixProtoproto) const = 0;
 #endif
 
+  // Write out the matrix as a sequence of (i,j,s) triplets. This
+  // format is useful for loading the matrix into MATLAB/octave as a
+  // sparse matrix.
+  virtual void ToTextFile(FILE* file) const = 0;
+
   // Accessors for the values array that stores the entries of the
   // sparse matrix. The exact interpreptation of the values of this
   // array depends on the particular kind of SparseMatrix being
index 7d7c3df9960dd7a8924e52af77554232259ee76b..247ab2e697b2cf691ec5b03b776293b29c1eef13 100644 (file)
@@ -295,5 +295,12 @@ TripletSparseMatrix* TripletSparseMatrix::CreateSparseDiagonalMatrix(
   return m;
 }
 
+void TripletSparseMatrix::ToTextFile(FILE* file) const {
+  CHECK_NOTNULL(file);
+  for (int i = 0; i < num_nonzeros_; ++i) {
+    fprintf(file, "% 10d % 10d %17f\n", rows_[i], cols_[i], values_[i]);
+  }
+}
+
 }  // namespace internal
 }  // namespace ceres
index 3c90a62fd20e1102db1552fdb9628e2cc7ea5063..300e74d0bbcb54f57c74b7a4ead282f5e2d840a9 100644 (file)
@@ -68,6 +68,7 @@ class TripletSparseMatrix : public SparseMatrix {
 #ifndef CERES_DONT_HAVE_PROTOCOL_BUFFERS
   virtual void ToProto(SparseMatrixProto *proto) const;
 #endif
+  virtual void ToTextFile(FILE* file) const;
   virtual int num_rows()        const { return num_rows_;     }
   virtual int num_cols()        const { return num_cols_;     }
   virtual int num_nonzeros()    const { return num_nonzeros_; }
diff --git a/extern/libmv/third_party/ceres/patches/msvc_isfinite.patch b/extern/libmv/third_party/ceres/patches/msvc_isfinite.patch
new file mode 100644 (file)
index 0000000..c3129d8
--- /dev/null
@@ -0,0 +1,15 @@
+diff --git a/internal/ceres/residual_block_utils.cc b/internal/ceres/residual_block_utils.cc
+index ed3499b..28e0313 100644
+--- a/internal/ceres/residual_block_utils.cc
++++ b/internal/ceres/residual_block_utils.cc
+@@ -40,6 +40,10 @@
+ #include "ceres/internal/eigen.h"
+ #include "ceres/internal/port.h"
++#ifdef _MSC_VER
++#  define isfinite _finite
++#endif
++
+ namespace ceres {
+ namespace internal {