← Back to branch summary

~ubuntu-branches/ubuntu/saucy/abyss/saucy

« back to all changes in this revision

Viewing changes to Graph/PopBubbles.h

  • Committer: Package Import Robot
  • Author(s): Shaun Jackman
  • Date: 2011-09-11 10:00:13 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: package-import@ubuntu.com-20110911100013-oa4m5fi036mjuwc8
Tags: 1.3.0-1
* New upstream release.
* Add libboost-graph-dev to Build-Depends.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Graph/PopBubbles.h
 
1
#ifndef POPBUBBLES_H
 
2
#define POPBUBBLES_H 1
 
3
 
 
4
#include "Common/Options.h"
 
5
#include "DepthFirstSearch.h"
 
6
#include "Graph.h"
 
7
#include <boost/unordered_map.hpp>
 
8
#include <algorithm>
 
9
#include <iostream>
 
10
#include <iterator>
 
11
#include <set>
 
12
#include <utility>
 
13
#include <vector>
 
14
 
 
15
/** Record a topological order of the vertices. */
 
16
template <typename OutIt>
 
17
struct TopoVisitor : public boost::default_dfs_visitor
 
18
{
 
19
        TopoVisitor(OutIt it) : m_it(it) { }
 
20
 
 
21
        template <typename Vertex, typename Graph>
 
22
        void finish_vertex(const Vertex& u, Graph&) { *m_it++ = u; }
 
23
 
 
24
  private:
 
25
        OutIt m_it;
 
26
};
 
27
 
 
28
/** Record a topological order of the vertices. */
 
29
template <typename Graph, typename It>
 
30
void topologicalSort(const Graph& g, It it)
 
31
{
 
32
        using boost::default_color_type;
 
33
        using boost::vector_property_map;
 
34
        typedef vector_property_map<
 
35
                default_color_type, ContigNodeIndexMap> ColorMap;
 
36
        depthFirstSearch(g, TopoVisitor<It>(it),
 
37
                        ColorMap(num_vertices(g)));
 
38
}
 
39
 
 
40
/** Return true if the specified sequence of vertices is a bubble. */
 
41
template <typename Graph, typename It>
 
42
bool isBubble(const Graph& g, It first, It last)
 
43
{
 
44
        typedef typename graph_traits<Graph>::adjacency_iterator Ait;
 
45
        typedef typename graph_traits<Graph>::vertex_descriptor V;
 
46
        assert(last - first > 1);
 
47
        if (last - first == 2)
 
48
                return false; // unambiguous edge
 
49
        if (*first == ~last[-1])
 
50
                return false; // palindrome
 
51
        std::set<V> targets(first, first + 1);
 
52
        for (It it = first; it != last - 1; ++it) {
 
53
                std::pair<Ait, Ait> adj = adjacent_vertices(*it, g);
 
54
                targets.insert(adj.first, adj.second);
 
55
        }
 
56
        std::set<V> sources(last - 1, last);
 
57
        for (It it = first + 1; it != last; ++it) {
 
58
                std::pair<Ait, Ait> adj = adjacent_vertices(~*it, g);
 
59
                transform(adj.first, adj.second,
 
60
                                inserter(sources, sources.end()),
 
61
                                std::mem_fun_ref(&V::operator~));
 
62
        }
 
63
        std::set<V> bubble(first, last);
 
64
        return sources == bubble && targets == bubble;
 
65
}
 
66
 
 
67
typedef std::vector<ContigNode> Bubble;
 
68
typedef std::vector<Bubble> Bubbles;
 
69
 
 
70
/** Discover bubbles. */
 
71
template <typename Graph>
 
72
Bubbles discoverBubbles(const Graph& g)
 
73
{
 
74
        typedef typename graph_traits<Graph>::vertex_descriptor V;
 
75
 
 
76
        std::vector<V> topo(num_vertices(g));
 
77
        topologicalSort(g, topo.rbegin());
 
78
 
 
79
        Bubbles bubbles;
 
80
        typedef typename std::vector<V>::const_iterator It;
 
81
        for (It first = topo.begin(); first != topo.end(); ++first) {
 
82
                int sum = out_degree(*first, g);
 
83
                if (sum < 2)
 
84
                        continue;
 
85
                if (opt::verbose > 3)
 
86
                        std::cerr << "* " << *first << '\n';
 
87
                for (It it = first + 1; it != topo.end(); ++it) {
 
88
                        unsigned indeg = in_degree(*it, g);
 
89
                        unsigned outdeg = out_degree(*it, g);
 
90
                        sum -= indeg;
 
91
 
 
92
                        if (opt::verbose > 3)
 
93
                                std::cerr << *it << '\t' << indeg << '\t' << outdeg
 
94
                                        << '\t' << sum
 
95
                                        << '\t' << sum + (int)outdeg << '\n';
 
96
 
 
97
                        if (indeg == 0 || sum < 0)
 
98
                                break;
 
99
                        if (sum == 0) {
 
100
                                It last = it + 1;
 
101
                                if (isBubble(g, first, last)) {
 
102
                                        if (opt::verbose > 3)
 
103
                                                std::cerr << "good\n";
 
104
                                        bubbles.push_back(Bubble(first, last));
 
105
                                        first = it - 1;
 
106
                                }
 
107
                                break;
 
108
                        }
 
109
 
 
110
                        if (outdeg == 0)
 
111
                                break;
 
112
                        sum += outdeg;
 
113
                }
 
114
        }
 
115
        return bubbles;
 
116
}
 
117
 
 
118
/** Return the length of the longest path through the bubble. */
 
119
template <typename Graph>
 
120
int longestPath(const Graph& g, const Bubble& topo)
 
121
{
 
122
        using boost::tie;
 
123
        typedef graph_traits<Graph> GTraits;
 
124
        typedef typename GTraits::edge_descriptor E;
 
125
        typedef typename GTraits::out_edge_iterator Eit;
 
126
        typedef typename GTraits::vertex_descriptor V;
 
127
 
 
128
        EdgeWeightMap<Graph> weight(g);
 
129
        /* GCC 4.4.6 has a bug that prevents ContigNode being used as the
 
130
         * key of a std::map. Possibly related to
 
131
         * http://gcc.gnu.org/bugzilla/show_bug.cgi?id=39390
 
132
         */
 
133
        boost::unordered_map<V, int> distance;
 
134
        distance[topo.front()] = 0;
 
135
        for (Bubble::const_iterator it = topo.begin();
 
136
                        it != topo.end(); ++it) {
 
137
                V u = *it;
 
138
                Eit eit, elast;
 
139
                for (tie(eit, elast) = out_edges(u, g); eit != elast; ++eit) {
 
140
                        E e = *eit;
 
141
                        V v = target(e, g);
 
142
                        distance[v] = std::max(
 
143
                                        distance[v], distance[u] + weight[e]);
 
144
                }
 
145
        }
 
146
        V v = topo.back();
 
147
        return distance[v] - g[v].length;
 
148
}
 
149
 
 
150
/** Scaffold over the bubble between vertices u and w.
 
151
 * Add an edge (u,w) with the distance property set to the length of
 
152
 * the largest branch of the bubble.
 
153
 */
 
154
template <typename Graph>
 
155
void scaffoldBubble(Graph& g, const Bubble& bubble)
 
156
{
 
157
        typedef graph_traits<Graph> GTraits;
 
158
        typedef typename GTraits::adjacency_iterator Ait;
 
159
        typedef typename GTraits::vertex_descriptor V;
 
160
        assert(bubble.size() > 2);
 
161
 
 
162
        V u = bubble.front(), w = bubble.back();
 
163
        if (edge(u, w, g).second) {
 
164
                // Already scaffolded.
 
165
                return;
 
166
        }
 
167
        assert(isBubble(g, bubble.begin(), bubble.end()));
 
168
 
 
169
        add_edge(u, w, std::max(longestPath(g, bubble), 1), g);
 
170
}
 
171
 
 
172
/** Replace each bubble in the graph with a single edge.
 
173
 * Remove the vertices in the bubbles from the graph.
 
174
 * @return the vertices that were removed from the graph
 
175
 */
 
176
template <typename Graph>
 
177
std::vector<typename graph_traits<Graph>::vertex_descriptor>
 
178
popBubbles(Graph& g)
 
179
{
 
180
        typedef typename graph_traits<Graph>::vertex_descriptor V;
 
181
        typedef std::vector<V> Vertices;
 
182
        Vertices popped;
 
183
        Bubbles bubbles = discoverBubbles(g);
 
184
        for (Bubbles::const_iterator it = bubbles.begin();
 
185
                        it != bubbles.end(); ++it) {
 
186
                scaffoldBubble(g, *it);
 
187
                popped.insert(popped.end(), it->begin() + 1, it->end() - 1);
 
188
        }
 
189
        for (typename Vertices::const_iterator it = popped.begin();
 
190
                        it != popped.end(); ++it) {
 
191
                V u = *it;
 
192
                clear_vertex(u, g);
 
193
                remove_vertex(u, g);
 
194
        }
 
195
        return popped;
 
196
}
 
197
 
 
198
#endif