Fix a bug in the reordering code.
When the user provides an ordering which starts at a non-zero group id,
or has gaps in the groups, then CAMD, the algorithm used to reorder
the program can crash or return garbage results.
The solution is to map the ordering into grouping constraints, and then
to re-number the groups to be contiguous using a call to
MapValuesToContiguousRange. This was already done for CAMD based
ordering for Schur type solvers, but was not done for SPARSE_NORMAL_CHOLESKY.
Thanks to Bernhard Zeisl for not only reporting the bug but also
providing a reproduction.
Change-Id: I5cfae222d701dfdb8e1bda7f0b4670a30417aa89
diff --git a/internal/ceres/parameter_block.h b/internal/ceres/parameter_block.h
index b50e116..cb7140d 100644
--- a/internal/ceres/parameter_block.h
+++ b/internal/ceres/parameter_block.h
@@ -238,9 +238,10 @@
}
std::string ToString() const {
- return StringPrintf("{ user_state=%p, state=%p, size=%d, "
+ return StringPrintf("{ this=%p, user_state=%p, state=%p, size=%d, "
"constant=%d, index=%d, state_offset=%d, "
"delta_offset=%d }",
+ this,
user_state_,
state_,
size_,
diff --git a/internal/ceres/reorder_program.cc b/internal/ceres/reorder_program.cc
index 90bd4ba..d0e8f32 100644
--- a/internal/ceres/reorder_program.cc
+++ b/internal/ceres/reorder_program.cc
@@ -132,6 +132,11 @@
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]);
ss.ConstrainedApproximateMinimumDegreeOrdering(block_jacobian_transpose,
&constraints[0],
ordering);
@@ -348,8 +353,8 @@
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,
+ // 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]);
@@ -535,6 +540,15 @@
const 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;
+ }
+
// Compute a block sparse presentation of J'.
scoped_ptr<TripletSparseMatrix> tsm_block_jacobian_transpose(
program->CreateJacobianBlockSparsityTranspose());
diff --git a/internal/ceres/reorder_program_test.cc b/internal/ceres/reorder_program_test.cc
index d1c07f8..5c0d516 100644
--- a/internal/ceres/reorder_program_test.cc
+++ b/internal/ceres/reorder_program_test.cc
@@ -36,6 +36,7 @@
#include "ceres/sized_cost_function.h"
#include "ceres/solver.h"
+#include "gmock/gmock.h"
#include "gtest/gtest.h"
namespace ceres {
@@ -168,5 +169,80 @@
EXPECT_EQ(parameter_blocks[2]->user_state(), &y);
}
+
+class ReorderProgramForSparseNormalCholeskyUsingSuiteSparseTest : public ::testing::Test {
+ protected:
+ void SetUp() {
+ problem_.AddResidualBlock(new UnaryCostFunction(), NULL, &x_);
+ problem_.AddResidualBlock(new BinaryCostFunction(), NULL, &z_, &x_);
+ problem_.AddResidualBlock(new BinaryCostFunction(), NULL, &z_, &y_);
+ problem_.AddResidualBlock(new UnaryCostFunction(), NULL, &z_);
+ problem_.AddResidualBlock(new BinaryCostFunction(), NULL, &x_, &y_);
+ problem_.AddResidualBlock(new UnaryCostFunction(), NULL, &y_);
+ };
+
+ void ComputeAndValidateOrdering(const ParameterBlockOrdering& linear_solver_ordering) {
+ Program* program = problem_.mutable_program();
+ vector<ParameterBlock*> unordered_parameter_blocks =
+ program->parameter_blocks();
+
+ std::string error;
+ EXPECT_TRUE(ReorderProgramForSparseNormalCholesky(
+ ceres::SUITE_SPARSE,
+ linear_solver_ordering,
+ program,
+ &error));
+ const vector<ParameterBlock*>& ordered_parameter_blocks =
+ program->parameter_blocks();
+ EXPECT_EQ(ordered_parameter_blocks.size(),
+ unordered_parameter_blocks.size());
+
+ EXPECT_THAT(unordered_parameter_blocks,
+ ::testing::UnorderedElementsAreArray(ordered_parameter_blocks));
+ }
+
+ ProblemImpl problem_;
+ double x_;
+ double y_;
+ double z_;
+};
+
+TEST_F(ReorderProgramForSparseNormalCholeskyUsingSuiteSparseTest, EverythingInGroupZero) {
+ ParameterBlockOrdering linear_solver_ordering;
+ linear_solver_ordering.AddElementToGroup(&x_, 0);
+ linear_solver_ordering.AddElementToGroup(&y_, 0);
+ linear_solver_ordering.AddElementToGroup(&z_, 0);
+
+ ComputeAndValidateOrdering(linear_solver_ordering);
+}
+
+TEST_F(ReorderProgramForSparseNormalCholeskyUsingSuiteSparseTest, ContiguousGroups) {
+ ParameterBlockOrdering linear_solver_ordering;
+ linear_solver_ordering.AddElementToGroup(&x_, 0);
+ linear_solver_ordering.AddElementToGroup(&y_, 1);
+ linear_solver_ordering.AddElementToGroup(&z_, 2);
+
+ ComputeAndValidateOrdering(linear_solver_ordering);
+}
+
+TEST_F(ReorderProgramForSparseNormalCholeskyUsingSuiteSparseTest, GroupsWithGaps) {
+ ParameterBlockOrdering linear_solver_ordering;
+ linear_solver_ordering.AddElementToGroup(&x_, 0);
+ linear_solver_ordering.AddElementToGroup(&y_, 2);
+ linear_solver_ordering.AddElementToGroup(&z_, 2);
+
+ ComputeAndValidateOrdering(linear_solver_ordering);
+}
+
+TEST_F(ReorderProgramForSparseNormalCholeskyUsingSuiteSparseTest,
+ NonContiguousStartingAtTwo) {
+ ParameterBlockOrdering linear_solver_ordering;
+ linear_solver_ordering.AddElementToGroup(&x_, 2);
+ linear_solver_ordering.AddElementToGroup(&y_, 4);
+ linear_solver_ordering.AddElementToGroup(&z_, 4);
+
+ ComputeAndValidateOrdering(linear_solver_ordering);
+}
+
} // namespace internal
} // namespace ceres
