Fix how Ceres calls CAMD.

CAMD requires that the id of the largest numbered elimination
group be less than the number of columns in the matrix.

This patch ensures that this is the case. Without this,
in certain cases its possible for CAMD to silently fail
while doing out of bounds access and then causing Ceres to fail.

Also add some logging about the problem size before and after
the reduced program has been created.

Change-Id: I0ea3c6572a7c29cbbf09afec9ba5b4f4d4b21a9b
diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc
index 06c65f9..83faa05 100644
--- a/internal/ceres/solver_impl.cc
+++ b/internal/ceres/solver_impl.cc
@@ -317,6 +317,16 @@
 void SolverImpl::Solve(const Solver::Options& options,
                        ProblemImpl* problem_impl,
                        Solver::Summary* summary) {
+  VLOG(2) << "Initial problem: "
+          << problem_impl->NumParameterBlocks()
+          << " parameter blocks, "
+          << problem_impl->NumParameters()
+          << " parameters,  "
+          << problem_impl->NumResidualBlocks()
+          << " residual blocks, "
+          << problem_impl->NumResiduals()
+          << " residuals.";
+
   if (options.minimizer_type == TRUST_REGION) {
     TrustRegionSolve(options, problem_impl, summary);
   } else {
@@ -1069,6 +1079,16 @@
     return NULL;
   }
 
+  VLOG(2) << "Reduced problem: "
+          << transformed_program->NumParameterBlocks()
+          << " parameter blocks, "
+          << transformed_program->NumParameters()
+          << " parameters,  "
+          << transformed_program->NumResidualBlocks()
+          << " residual blocks, "
+          << transformed_program->NumResiduals()
+          << " residuals.";
+
   if (transformed_program->NumParameterBlocks() == 0) {
     LOG(WARNING) << "No varying parameter blocks to optimize; "
                  << "bailing early.";
@@ -1617,6 +1637,11 @@
               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].
+    SolverImpl::CompactifyArray(&constraints);
+
     // Set the offsets and index for CreateJacobianSparsityTranspose.
     program->SetParameterOffsetsAndIndex();
     // Compute a block sparse presentation of J'.
@@ -1740,5 +1765,20 @@
   return true;
 }
 
+void SolverImpl::CompactifyArray(vector<int>* array_ptr) {
+  vector<int>& array = *array_ptr;
+  const set<int> unique_group_ids(array.begin(), array.end());
+  map<int, int> group_id_map;
+  for (set<int>::const_iterator it = unique_group_ids.begin();
+       it != unique_group_ids.end();
+       ++it) {
+    InsertOrDie(&group_id_map, *it, group_id_map.size());
+  }
+
+  for (int i = 0; i < array.size(); ++i) {
+    array[i] = group_id_map[array[i]];
+  }
+}
+
 }  // namespace internal
 }  // namespace ceres
diff --git a/internal/ceres/solver_impl.h b/internal/ceres/solver_impl.h
index ebfb813..2b7ca3e 100644
--- a/internal/ceres/solver_impl.h
+++ b/internal/ceres/solver_impl.h
@@ -208,6 +208,13 @@
       ParameterBlockOrdering* parameter_block_ordering,
       Program* program,
       string* error);
+
+  // 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
diff --git a/internal/ceres/solver_impl_test.cc b/internal/ceres/solver_impl_test.cc
index 4f1e549..2958e84 100644
--- a/internal/ceres/solver_impl_test.cc
+++ b/internal/ceres/solver_impl_test.cc
@@ -991,5 +991,52 @@
   EXPECT_EQ((expected_dense_jacobian - actual_dense_jacobian).norm(), 0.0);
 }
 
+TEST(CompactifyArray, ContiguousEntries) {
+  vector<int> array;
+  array.push_back(0);
+  array.push_back(1);
+  vector<int> expected = array;
+  SolverImpl::CompactifyArray(&array);
+  EXPECT_EQ(array, expected);
+  array.clear();
+
+  array.push_back(1);
+  array.push_back(0);
+  expected = array;
+  SolverImpl::CompactifyArray(&array);
+  EXPECT_EQ(array, expected);
+}
+
+TEST(CompactifyArray, NonContiguousEntries) {
+  vector<int> array;
+  array.push_back(0);
+  array.push_back(2);
+  vector<int> expected;
+  expected.push_back(0);
+  expected.push_back(1);
+  SolverImpl::CompactifyArray(&array);
+  EXPECT_EQ(array, expected);
+}
+
+TEST(CompactifyArray, ) {
+  vector<int> array;
+  array.push_back(3);
+  array.push_back(1);
+  array.push_back(0);
+  array.push_back(0);
+  array.push_back(0);
+  array.push_back(5);
+  vector<int> expected;
+  expected.push_back(2);
+  expected.push_back(1);
+  expected.push_back(0);
+  expected.push_back(0);
+  expected.push_back(0);
+  expected.push_back(3);
+
+  SolverImpl::CompactifyArray(&array);
+  EXPECT_EQ(array, expected);
+}
+
 }  // namespace internal
 }  // namespace ceres