1// Ceres Solver - A fast non-linear least squares minimizer 2// Copyright 2010, 2011, 2012 Google Inc. All rights reserved. 3// http://code.google.com/p/ceres-solver/ 4// 5// Redistribution and use in source and binary forms, with or without 6// modification, are permitted provided that the following conditions are met: 7// 8// * Redistributions of source code must retain the above copyright notice, 9// this list of conditions and the following disclaimer. 10// * Redistributions in binary form must reproduce the above copyright notice, 11// this list of conditions and the following disclaimer in the documentation 12// and/or other materials provided with the distribution. 13// * Neither the name of Google Inc. nor the names of its contributors may be 14// used to endorse or promote products derived from this software without 15// specific prior written permission. 16// 17// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 18// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 19// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 20// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 21// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 22// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 23// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 24// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 25// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 26// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 27// POSSIBILITY OF SUCH DAMAGE. 28// 29// Author: sameeragarwal@google.com (Sameer Agarwal) 30// 31// Various algorithms that operate on undirected graphs. 32 33#ifndef CERES_INTERNAL_GRAPH_ALGORITHMS_H_ 34#define CERES_INTERNAL_GRAPH_ALGORITHMS_H_ 35 36#include <vector> 37#include <glog/logging.h> 38#include "ceres/collections_port.h" 39#include "ceres/graph.h" 40 41namespace ceres { 42namespace internal { 43 44// Compare two vertices of a graph by their degrees. 45template <typename Vertex> 46class VertexDegreeLessThan { 47 public: 48 explicit VertexDegreeLessThan(const Graph<Vertex>& graph) 49 : graph_(graph) {} 50 51 bool operator()(const Vertex& lhs, const Vertex& rhs) const { 52 if (graph_.Neighbors(lhs).size() == graph_.Neighbors(rhs).size()) { 53 return lhs < rhs; 54 } 55 return graph_.Neighbors(lhs).size() < graph_.Neighbors(rhs).size(); 56 } 57 58 private: 59 const Graph<Vertex>& graph_; 60}; 61 62// Order the vertices of a graph using its (approximately) largest 63// independent set, where an independent set of a graph is a set of 64// vertices that have no edges connecting them. The maximum 65// independent set problem is NP-Hard, but there are effective 66// approximation algorithms available. The implementation here uses a 67// breadth first search that explores the vertices in order of 68// increasing degree. The same idea is used by Saad & Li in "MIQR: A 69// multilevel incomplete QR preconditioner for large sparse 70// least-squares problems", SIMAX, 2007. 71// 72// Given a undirected graph G(V,E), the algorithm is a greedy BFS 73// search where the vertices are explored in increasing order of their 74// degree. The output vector ordering contains elements of S in 75// increasing order of their degree, followed by elements of V - S in 76// increasing order of degree. The return value of the function is the 77// cardinality of S. 78template <typename Vertex> 79int IndependentSetOrdering(const Graph<Vertex>& graph, 80 vector<Vertex>* ordering) { 81 const HashSet<Vertex>& vertices = graph.vertices(); 82 const int num_vertices = vertices.size(); 83 84 CHECK_NOTNULL(ordering); 85 ordering->clear(); 86 ordering->reserve(num_vertices); 87 88 // Colors for labeling the graph during the BFS. 89 const char kWhite = 0; 90 const char kGrey = 1; 91 const char kBlack = 2; 92 93 // Mark all vertices white. 94 HashMap<Vertex, char> vertex_color; 95 vector<Vertex> vertex_queue; 96 for (typename HashSet<Vertex>::const_iterator it = vertices.begin(); 97 it != vertices.end(); 98 ++it) { 99 vertex_color[*it] = kWhite; 100 vertex_queue.push_back(*it); 101 } 102 103 104 sort(vertex_queue.begin(), vertex_queue.end(), 105 VertexDegreeLessThan<Vertex>(graph)); 106 107 // Iterate over vertex_queue. Pick the first white vertex, add it 108 // to the independent set. Mark it black and its neighbors grey. 109 for (int i = 0; i < vertex_queue.size(); ++i) { 110 const Vertex& vertex = vertex_queue[i]; 111 if (vertex_color[vertex] != kWhite) { 112 continue; 113 } 114 115 ordering->push_back(vertex); 116 vertex_color[vertex] = kBlack; 117 const HashSet<Vertex>& neighbors = graph.Neighbors(vertex); 118 for (typename HashSet<Vertex>::const_iterator it = neighbors.begin(); 119 it != neighbors.end(); 120 ++it) { 121 vertex_color[*it] = kGrey; 122 } 123 } 124 125 int independent_set_size = ordering->size(); 126 127 // Iterate over the vertices and add all the grey vertices to the 128 // ordering. At this stage there should only be black or grey 129 // vertices in the graph. 130 for (typename vector<Vertex>::const_iterator it = vertex_queue.begin(); 131 it != vertex_queue.end(); 132 ++it) { 133 const Vertex vertex = *it; 134 DCHECK(vertex_color[vertex] != kWhite); 135 if (vertex_color[vertex] != kBlack) { 136 ordering->push_back(vertex); 137 } 138 } 139 140 CHECK_EQ(ordering->size(), num_vertices); 141 return independent_set_size; 142} 143 144// Find the connected component for a vertex implemented using the 145// find and update operation for disjoint-set. Recursively traverse 146// the disjoint set structure till you reach a vertex whose connected 147// component has the same id as the vertex itself. Along the way 148// update the connected components of all the vertices. This updating 149// is what gives this data structure its efficiency. 150template <typename Vertex> 151Vertex FindConnectedComponent(const Vertex& vertex, 152 HashMap<Vertex, Vertex>* union_find) { 153 typename HashMap<Vertex, Vertex>::iterator it = union_find->find(vertex); 154 DCHECK(it != union_find->end()); 155 if (it->second != vertex) { 156 it->second = FindConnectedComponent(it->second, union_find); 157 } 158 159 return it->second; 160} 161 162// Compute a degree two constrained Maximum Spanning Tree/forest of 163// the input graph. Caller owns the result. 164// 165// Finding degree 2 spanning tree of a graph is not always 166// possible. For example a star graph, i.e. a graph with n-nodes 167// where one node is connected to the other n-1 nodes does not have 168// a any spanning trees of degree less than n-1.Even if such a tree 169// exists, finding such a tree is NP-Hard. 170 171// We get around both of these problems by using a greedy, degree 172// constrained variant of Kruskal's algorithm. We start with a graph 173// G_T with the same vertex set V as the input graph G(V,E) but an 174// empty edge set. We then iterate over the edges of G in decreasing 175// order of weight, adding them to G_T if doing so does not create a 176// cycle in G_T} and the degree of all the vertices in G_T remains 177// bounded by two. This O(|E|) algorithm results in a degree-2 178// spanning forest, or a collection of linear paths that span the 179// graph G. 180template <typename Vertex> 181Graph<Vertex>* 182Degree2MaximumSpanningForest(const Graph<Vertex>& graph) { 183 // Array of edges sorted in decreasing order of their weights. 184 vector<pair<double, pair<Vertex, Vertex> > > weighted_edges; 185 Graph<Vertex>* forest = new Graph<Vertex>(); 186 187 // Disjoint-set to keep track of the connected components in the 188 // maximum spanning tree. 189 HashMap<Vertex, Vertex> disjoint_set; 190 191 // Sort of the edges in the graph in decreasing order of their 192 // weight. Also add the vertices of the graph to the Maximum 193 // Spanning Tree graph and set each vertex to be its own connected 194 // component in the disjoint_set structure. 195 const HashSet<Vertex>& vertices = graph.vertices(); 196 for (typename HashSet<Vertex>::const_iterator it = vertices.begin(); 197 it != vertices.end(); 198 ++it) { 199 const Vertex vertex1 = *it; 200 forest->AddVertex(vertex1, graph.VertexWeight(vertex1)); 201 disjoint_set[vertex1] = vertex1; 202 203 const HashSet<Vertex>& neighbors = graph.Neighbors(vertex1); 204 for (typename HashSet<Vertex>::const_iterator it2 = neighbors.begin(); 205 it2 != neighbors.end(); 206 ++it2) { 207 const Vertex vertex2 = *it2; 208 if (vertex1 >= vertex2) { 209 continue; 210 } 211 const double weight = graph.EdgeWeight(vertex1, vertex2); 212 weighted_edges.push_back(make_pair(weight, make_pair(vertex1, vertex2))); 213 } 214 } 215 216 // The elements of this vector, are pairs<edge_weight, 217 // edge>. Sorting it using the reverse iterators gives us the edges 218 // in decreasing order of edges. 219 sort(weighted_edges.rbegin(), weighted_edges.rend()); 220 221 // Greedily add edges to the spanning tree/forest as long as they do 222 // not violate the degree/cycle constraint. 223 for (int i =0; i < weighted_edges.size(); ++i) { 224 const pair<Vertex, Vertex>& edge = weighted_edges[i].second; 225 const Vertex vertex1 = edge.first; 226 const Vertex vertex2 = edge.second; 227 228 // Check if either of the vertices are of degree 2 already, in 229 // which case adding this edge will violate the degree 2 230 // constraint. 231 if ((forest->Neighbors(vertex1).size() == 2) || 232 (forest->Neighbors(vertex2).size() == 2)) { 233 continue; 234 } 235 236 // Find the id of the connected component to which the two 237 // vertices belong to. If the id is the same, it means that the 238 // two of them are already connected to each other via some other 239 // vertex, and adding this edge will create a cycle. 240 Vertex root1 = FindConnectedComponent(vertex1, &disjoint_set); 241 Vertex root2 = FindConnectedComponent(vertex2, &disjoint_set); 242 243 if (root1 == root2) { 244 continue; 245 } 246 247 // This edge can be added, add an edge in either direction with 248 // the same weight as the original graph. 249 const double edge_weight = graph.EdgeWeight(vertex1, vertex2); 250 forest->AddEdge(vertex1, vertex2, edge_weight); 251 forest->AddEdge(vertex2, vertex1, edge_weight); 252 253 // Connected the two connected components by updating the 254 // disjoint_set structure. Always connect the connected component 255 // with the greater index with the connected component with the 256 // smaller index. This should ensure shallower trees, for quicker 257 // lookup. 258 if (root2 < root1) { 259 std::swap(root1, root2); 260 }; 261 262 disjoint_set[root2] = root1; 263 } 264 return forest; 265} 266 267} // namespace internal 268} // namespace ceres 269 270#endif // CERES_INTERNAL_GRAPH_ALGORITHMS_H_ 271