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/
5
// Redistribution and use in source and binary forms, with or without
6
// modification, are permitted provided that the following conditions are met:
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.
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.
29
// Author: sameeragarwal@google.com (Sameer Agarwal)
31
// Various algorithms that operate on undirected graphs.
33
#ifndef CERES_INTERNAL_GRAPH_ALGORITHMS_H_
34
#define CERES_INTERNAL_GRAPH_ALGORITHMS_H_
37
#include <glog/logging.h>
38
#include "ceres/collections_port.h"
39
#include "ceres/graph.h"
44
// Compare two vertices of a graph by their degrees.
45
template <typename Vertex>
46
class VertexDegreeLessThan {
48
explicit VertexDegreeLessThan(const Graph<Vertex>& graph)
51
bool operator()(const Vertex& lhs, const Vertex& rhs) const {
52
if (graph_.Neighbors(lhs).size() == graph_.Neighbors(rhs).size()) {
55
return graph_.Neighbors(lhs).size() < graph_.Neighbors(rhs).size();
59
const Graph<Vertex>& graph_;
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.
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
78
template <typename Vertex>
79
int IndependentSetOrdering(const Graph<Vertex>& graph,
80
vector<Vertex>* ordering) {
81
const HashSet<Vertex>& vertices = graph.vertices();
82
const int num_vertices = vertices.size();
84
CHECK_NOTNULL(ordering);
86
ordering->reserve(num_vertices);
88
// Colors for labeling the graph during the BFS.
89
const char kWhite = 0;
91
const char kBlack = 2;
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();
99
vertex_color[*it] = kWhite;
100
vertex_queue.push_back(*it);
104
sort(vertex_queue.begin(), vertex_queue.end(),
105
VertexDegreeLessThan<Vertex>(graph));
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) {
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();
121
vertex_color[*it] = kGrey;
125
int independent_set_size = ordering->size();
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();
133
const Vertex vertex = *it;
134
DCHECK(vertex_color[vertex] != kWhite);
135
if (vertex_color[vertex] != kBlack) {
136
ordering->push_back(vertex);
140
CHECK_EQ(ordering->size(), num_vertices);
141
return independent_set_size;
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.
150
template <typename Vertex>
151
Vertex 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);
162
// Compute a degree two constrained Maximum Spanning Tree/forest of
163
// the input graph. Caller owns the result.
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.
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
180
template <typename Vertex>
182
Degree2MaximumSpanningForest(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>();
187
// Disjoint-set to keep track of the connected components in the
188
// maximum spanning tree.
189
HashMap<Vertex, Vertex> disjoint_set;
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();
199
const Vertex vertex1 = *it;
200
forest->AddVertex(vertex1, graph.VertexWeight(vertex1));
201
disjoint_set[vertex1] = vertex1;
203
const HashSet<Vertex>& neighbors = graph.Neighbors(vertex1);
204
for (typename HashSet<Vertex>::const_iterator it2 = neighbors.begin();
205
it2 != neighbors.end();
207
const Vertex vertex2 = *it2;
208
if (vertex1 >= vertex2) {
211
const double weight = graph.EdgeWeight(vertex1, vertex2);
212
weighted_edges.push_back(make_pair(weight, make_pair(vertex1, vertex2)));
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());
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;
228
// Check if either of the vertices are of degree 2 already, in
229
// which case adding this edge will violate the degree 2
231
if ((forest->Neighbors(vertex1).size() == 2) ||
232
(forest->Neighbors(vertex2).size() == 2)) {
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);
243
if (root1 == root2) {
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);
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
259
std::swap(root1, root2);
262
disjoint_set[root2] = root1;
267
} // namespace internal
270
#endif // CERES_INTERNAL_GRAPH_ALGORITHMS_H_
added
removed
internal/ceres/graph_algorithms.h
