graph_distance.cc 16.6 KB
Newer Older
Tiago Peixoto's avatar
Tiago Peixoto committed
1
2
3
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2015 Tiago de Paula Peixoto <tiago@skewed.de>
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 3
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.

#include "graph.hh"
#include "graph_filtering.hh"
#include "graph_properties.hh"
#include "graph_selectors.hh"
Tiago Peixoto's avatar
Tiago Peixoto committed
22
#include "graph_python_interface.hh"
23
#include "numpy_bind.hh"
24
#include "hash_map_wrap.hh"
25
26

#include <boost/graph/breadth_first_search.hpp>
27
28
#include <boost/graph/dijkstra_shortest_paths_no_color_map.hpp>
#include <boost/graph/bellman_ford_shortest_paths.hpp>
29
#include <boost/python/stl_iterator.hpp>
30
31
#include <boost/python.hpp>

Tiago Peixoto's avatar
Tiago Peixoto committed
32
33
34
35
#ifdef HAVE_BOOST_COROUTINE
#include <boost/coroutine/all.hpp>
#endif // HAVE_BOOST_COROUTINE

36
37
38
39
40
41
42
43
44
45
46
using namespace std;
using namespace boost;
using namespace graph_tool;

struct stop_search {};

template <class DistMap, class PredMap>
class bfs_max_visitor:
    public boost::bfs_visitor<null_visitor>
{
public:
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
    bfs_max_visitor(DistMap dist_map, PredMap pred, size_t max_dist,
                    size_t source, size_t target)
        : _dist_map(dist_map), _pred(pred), _max_dist(max_dist),
          _source(source), _target(target), _dist(0) {}

    template <class Graph>
    void initialize_vertex(typename graph_traits<Graph>::vertex_descriptor v,
                           Graph&)
    {
        typedef typename property_traits<DistMap>::value_type dist_t;
        dist_t inf = std::is_floating_point<dist_t>::value ?
            numeric_limits<dist_t>::infinity() :
            numeric_limits<dist_t>::max();
        _dist_map[v] = (v == _source) ? 0 : inf;
        _pred[v] = v;
    }
63
64

    template <class Graph>
65
66
    void tree_edge(typename graph_traits<Graph>::edge_descriptor e,
                   Graph& g)
67
68
69
70
    {
        _pred[target(e,g)] = source(e,g);
    }

71
72
73
74
    template <class Graph>
    void examine_vertex(typename graph_traits<Graph>::vertex_descriptor v,
                        Graph&)
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
75
76
        typedef typename property_traits<DistMap>::value_type val_t;
        if ( _dist_map[v] > val_t(_max_dist))
77
78
79
            throw stop_search();
    }

80
81
    template <class Graph>
    void discover_vertex(typename graph_traits<Graph>::vertex_descriptor v,
82
                         Graph&)
83
    {
84
        if (size_t(_pred[v]) == v)
85
            return;
86
        _dist_map[v] = _dist_map[_pred[v]] + 1;
87
88
        if (v == _target)
            throw stop_search();
89
90
91
92
93
94
    }

private:
    DistMap _dist_map;
    PredMap _pred;
    size_t _max_dist;
95
    size_t _source;
96
    size_t _target;
97
98
99
    size_t _dist;
};

100
101
102
103
104
template <class DistMap, class PredMap>
class bfs_max_multiple_targets_visitor:
    public boost::bfs_visitor<null_visitor>
{
public:
105
    bfs_max_multiple_targets_visitor(DistMap dist_map, PredMap pred,
106
                                     size_t max_dist, size_t source,
107
                                     gt_hash_set<std::size_t> target)
108
109
110
111
112
113
114
115
116
117
118
119
120
121
        : _dist_map(dist_map), _pred(pred), _max_dist(max_dist),
          _source(source), _target(target), _dist(0) {}

    template <class Graph>
    void initialize_vertex(typename graph_traits<Graph>::vertex_descriptor v,
                           Graph&)
    {
        typedef typename property_traits<DistMap>::value_type dist_t;
        dist_t inf = std::is_floating_point<dist_t>::value ?
            numeric_limits<dist_t>::infinity() :
            numeric_limits<dist_t>::max();
        _dist_map[v] = (v == _source) ? 0 : inf;
        _pred[v] = v;
    }
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146

    template <class Graph>
    void tree_edge(typename graph_traits<Graph>::edge_descriptor e,
                   Graph& g)
    {
        _pred[target(e,g)] = source(e,g);
    }

    template <class Graph>
    void examine_vertex(typename graph_traits<Graph>::vertex_descriptor v,
                        Graph&)
    {
        typedef typename property_traits<DistMap>::value_type val_t;
        if ( _dist_map[v] > val_t(_max_dist))
            throw stop_search();
    }

    template <class Graph>
    void discover_vertex(typename graph_traits<Graph>::vertex_descriptor v,
                         Graph&)
    {
        if (size_t(_pred[v]) == v)
            return;
        _dist_map[v] = _dist_map[_pred[v]] + 1;

147
148
        auto iter = _target.find(v);
        if (iter != _target.end())
149
        {
150
            _target.erase(iter);
151
152
153
154
155
156
157
158
159
            if (_target.empty())
                throw stop_search();
        };
    }

private:
    DistMap _dist_map;
    PredMap _pred;
    size_t _max_dist;
160
    size_t _source;
161
    gt_hash_set<std::size_t> _target;
162
163
164
165
    size_t _dist;
};


166
template <class DistMap>
167
168
169
170
class djk_max_visitor:
    public boost::dijkstra_visitor<null_visitor>
{
public:
171
    djk_max_visitor(DistMap dist_map,
172
173
                    typename property_traits<DistMap>::value_type max_dist,
                    size_t target)
174
        : _dist_map(dist_map), _max_dist(max_dist), _target(target) {}
175

176

177
    template <class Graph>
178
179
    void examine_vertex(typename graph_traits<Graph>::vertex_descriptor u,
                        Graph&)
180
    {
181
        if (_dist_map[u] > _max_dist)
182
            throw stop_search();
183

184
185
        if (u == _target)
            throw stop_search();
186
187
    }

188

189
190
191
private:
    DistMap _dist_map;
    typename property_traits<DistMap>::value_type _max_dist;
192
    size_t _target;
193
194
195
};


196
197
198
199
200
201
202
template <class DistMap>
class djk_max_multiple_targets_visitor:
    public boost::dijkstra_visitor<null_visitor>
{
public:
    djk_max_multiple_targets_visitor(DistMap dist_map,
                                     typename property_traits<DistMap>::value_type max_dist, 
203
                                     gt_hash_set<std::size_t> target)
204
205
206
207
208
209
210
211
212
        : _dist_map(dist_map), _max_dist(max_dist), _target(target) {}


    template <class Graph>
    void examine_vertex(typename graph_traits<Graph>::vertex_descriptor u,
                        Graph&)
    {
        if (_dist_map[u] > _max_dist)
            throw stop_search();
213

214
215
        auto iter = _target.find(u);
        if (iter != _target.end())
216
        {
217
            _target.erase(iter);
218
219
220
221
222
223
224
225
226
            if (_target.empty())
                throw stop_search();
        };
    }


private:
    DistMap _dist_map;
    typename property_traits<DistMap>::value_type _max_dist;
227
    gt_hash_set<std::size_t> _target;
228
229
};

230
231
struct do_bfs_search
{
232
    template <class Graph, class VertexIndexMap, class DistMap, class PredMap>
233
234
    void operator()(const Graph& g, size_t source,
                    boost::python::object otarget_list,
235
236
                    VertexIndexMap vertex_index, DistMap dist_map,
                    PredMap pred_map, long double max_dist) const
237
238
    {
        typedef typename property_traits<DistMap>::value_type dist_t;
239

240
        auto target_list = get_array<int64_t, 1>(otarget_list);
241
242
        gt_hash_set<std::size_t> tgt(target_list.begin(),
                                     target_list.end());
243

244
245
246
        dist_t inf = std::is_floating_point<dist_t>::value ?
            numeric_limits<dist_t>::infinity() :
            numeric_limits<dist_t>::max();
247

248
        dist_t max_d = (max_dist > 0) ? max_dist : inf;
249

250
        unchecked_vector_property_map<boost::default_color_type, VertexIndexMap>
251
        color_map(vertex_index, num_vertices(g));
252
        try
253
        {
254
255
            if (tgt.size() <= 1)
            {
256
257
258
                size_t target = tgt.empty() ?
                    graph_traits<GraphInterface::multigraph_t>::null_vertex() :
                    *tgt.begin();
259
260
                breadth_first_search(g, vertex(source, g),
                                     visitor(bfs_max_visitor<DistMap, PredMap>
261
262
                                             (dist_map, pred_map, max_d,
                                              source, target)).
263
264
265
266
267
268
269
                                     vertex_index_map(vertex_index).
                                     color_map(color_map));
            }
            else
            {
                breadth_first_search(g, vertex(source, g),
                                     visitor(bfs_max_multiple_targets_visitor<DistMap, PredMap>
270
271
                                             (dist_map, pred_map, max_d,
                                              source, tgt)).
272
273
274
275
                                     vertex_index_map(vertex_index).
                                     color_map(color_map));
            }

276
277
278
279
280
281
282
        }
        catch (stop_search&) {}
    }
};

struct do_djk_search
{
283
284
    template <class Graph, class VertexIndexMap, class DistMap, class PredMap,
              class WeightMap>
285
286
    void operator()(const Graph& g, size_t source,
                    boost::python::object otarget_list,
287
288
                    VertexIndexMap vertex_index, DistMap dist_map,
                    PredMap pred_map, WeightMap weight, long double max_dist) const
289
    {
290
        auto target_list = get_array<int64_t, 1>(otarget_list);
291
292
        typedef typename property_traits<DistMap>::value_type dist_t;
        dist_t max_d = (max_dist > 0) ?
293
294
295
            max_dist : (std::is_floating_point<dist_t>::value ?
                        numeric_limits<dist_t>::infinity() :
                        numeric_limits<dist_t>::max());
296

297
298
        gt_hash_set<std::size_t> tgt(target_list.begin(),
                                     target_list.end());
299

300
301
302
303
        dist_t inf = (std::is_floating_point<dist_t>::value) ?
            numeric_limits<dist_t>::infinity() :
            numeric_limits<dist_t>::max();

304
        int i, N = num_vertices(g);
305
306
        #pragma omp parallel for default(shared) private(i) \
            schedule(runtime) if (N > 100)
307
        for (i = 0; i < N; ++i)
308
            dist_map[i] = inf;
309
310
        dist_map[source] = 0;

311
        try
312
        {
313
314
            if (tgt.size() <= 1)
            {
315
316
317
                size_t target = tgt.empty() ?
                    graph_traits<GraphInterface::multigraph_t>::null_vertex() :
                    *tgt.begin();
318
319
320
321
322
323
324
325
326
                dijkstra_shortest_paths_no_color_map
                    (g, vertex(source, g),
                     weight_map(weight).
                     distance_map(dist_map).
                     vertex_index_map(vertex_index).
                     predecessor_map(pred_map).
                     distance_inf(inf).
                     visitor(djk_max_visitor<DistMap>
                             (dist_map, max_d, target)));
327
328
329
            }
            else
            {
330
331
332
333
334
335
336
337
338
                dijkstra_shortest_paths_no_color_map
                    (g, vertex(source, g),
                     weight_map(weight).
                     distance_map(dist_map).
                     vertex_index_map(vertex_index).
                     predecessor_map(pred_map).
                     distance_inf(inf).
                     visitor(djk_max_multiple_targets_visitor<DistMap>
                             (dist_map, max_d, tgt)));
339
340
            }

341
342
343
344
345
        }
        catch (stop_search&) {}
    }
};

346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
struct do_bf_search
{
    template <class Graph, class DistMap, class PredMap, class WeightMap>
    void operator()(const Graph& g, size_t source, DistMap dist_map,
                    PredMap pred_map, WeightMap weight) const
    {
        bool ret = bellman_ford_shortest_paths(g, root_vertex(source).
                                               predecessor_map(pred_map).
                                               distance_map(dist_map).
                                               weight_map(weight));
        if (!ret)
            throw ValueException("Graph contains negative loops");

        // consistency with dijkstra
        typedef typename property_traits<DistMap>::value_type dist_t;
        if (std::is_floating_point<dist_t>::value)
        {
            for (auto v : vertices_range(g))
            {
                if (dist_map[v] == numeric_limits<dist_t>::max())
                    dist_map[v] = numeric_limits<dist_t>::infinity();
            }
        }
    }
};

372
373
void get_dists(GraphInterface& gi, size_t source, boost::python::object tgt,
               boost::any dist_map, boost::any weight, boost::any pred_map,
374
               long double max_dist, bool bf)
375
{
376
377
378
379
380
    typedef property_map_type
        ::apply<int64_t, GraphInterface::vertex_index_map_t>::type pred_map_t;

    pred_map_t pmap = any_cast<pred_map_t>(pred_map);

381
382
383
    if (weight.empty())
    {
        run_action<>()
384
385
            (gi, std::bind(do_bfs_search(), placeholders::_1, source, tgt, gi.get_vertex_index(),
                           placeholders::_2, pmap.get_unchecked(num_vertices(gi.get_graph())),
Tiago Peixoto's avatar
Tiago Peixoto committed
386
                           max_dist),
387
             writable_vertex_scalar_properties())
388
389
390
391
            (dist_map);
    }
    else
    {
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
        if (bf)
        {
            run_action<>()
                (gi, std::bind(do_bf_search(), placeholders::_1, source,
                               placeholders::_2, pmap.get_unchecked(num_vertices(gi.get_graph())),
                               placeholders::_3),
                 writable_vertex_scalar_properties(),
                 edge_scalar_properties())
                (dist_map, weight);
        }
        else
        {
            run_action<>()
                (gi, std::bind(do_djk_search(), placeholders::_1, source, tgt, gi.get_vertex_index(),
                               placeholders::_2, pmap.get_unchecked(num_vertices(gi.get_graph())),
                               placeholders::_3, max_dist),
                 writable_vertex_scalar_properties(),
                 edge_scalar_properties())
                (dist_map, weight);
        }
412
413
414
    }
}

Tiago Peixoto's avatar
Tiago Peixoto committed
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
template <class Graph, class Dist, class Pred, class Preds>
void get_all_preds(Graph g, Dist dist, Pred pred, Preds preds)
{
    int i, N = num_vertices(g);
    #pragma omp parallel for default(shared) private(i)     \
        schedule(runtime) if (N > 100)
    for (i = 0; i < N; ++i)
    {
        size_t v = vertex(i, g);
        if (v == graph_traits<Graph>::null_vertex())
            continue;
        if (size_t(pred[v]) == v)
            continue;
        auto d = dist[pred[v]];
        for (auto e : in_or_out_edges_range(v, g))
        {
            auto u = boost::is_directed(g) ? source(e, g) : target(e, g);
            if (dist[u] == d)
                preds[v].push_back(u);
        }
    }
};

void do_get_all_preds(GraphInterface& gi, boost::any adist,
                      boost::any apred, boost::any apreds)
{
    typedef property_map_type
        ::apply<int64_t, GraphInterface::vertex_index_map_t>::type pred_map_t;
    typedef property_map_type
        ::apply<vector<int64_t>, GraphInterface::vertex_index_map_t>::type preds_map_t;

    pred_map_t pred = any_cast<pred_map_t>(apred);
    preds_map_t preds = any_cast<preds_map_t>(apreds);

    run_action<>()
        (gi, [&](auto& g, auto dist)
             {get_all_preds(g, dist, pred.get_unchecked(num_vertices(g)),
                            preds.get_unchecked(num_vertices(g)));},
         vertex_scalar_properties())(adist);
}


template <class Pred, class Yield>
void get_all_paths(size_t s, size_t t, Pred pred, Yield& yield)
{
    vector<size_t> path;
    vector<pair<size_t, size_t>> stack = {{t, 0}};
    while (!stack.empty())
    {
        size_t v, i;
        std::tie(v, i) = stack.back();
        if (v == s)
        {
            path.clear();
            for (auto iter = stack.rbegin(); iter != stack.rend(); ++iter)
                path.push_back(iter->first);
            yield(wrap_vector_owned<size_t>(path));
        }
        if (pred[v].size() > i)
        {
            stack.emplace_back(pred[v][i], 0);
        }
        else
        {
            stack.pop_back();
            if (!stack.empty())
                ++stack.back().second;
        }
    }
};

python::object do_get_all_paths(GraphInterface& gi, size_t s, size_t t,
                                boost::any apred)
{
#ifdef HAVE_BOOST_COROUTINE
    auto dispatch = [&](auto& yield)
        {
            run_action<>()
                (gi, [&](auto&, auto pred) {get_all_paths(s, t, pred, yield);},
                 vertex_scalar_vector_properties())(apred);
        };
    return python::object(CoroGenerator(dispatch));
#else
    throw GraphException("This functionality is not available because boost::coroutine was not found at compile-time");
#endif // HAVE_BOOST_COROUTINE
}

502
503
504
void export_dists()
{
    python::def("get_dists", &get_dists);
Tiago Peixoto's avatar
Tiago Peixoto committed
505
506
    python::def("get_all_preds", &do_get_all_preds);
    python::def("get_all_paths", &do_get_all_paths);
507
};