graph_blockmodel.cc 30.8 KB
Newer Older
1
2
// graph-tool -- a general graph modification and manipulation thingy
//
Tiago Peixoto's avatar
Tiago Peixoto committed
3
// 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
//
// 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/>.

18
#define BOOST_PYTHON_MAX_ARITY 40
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#include <boost/python.hpp>
#include <cmath>
#include <iostream>

#include "numpy_bind.hh"

#include <boost/python.hpp>
#include <boost/python/suite/indexing/vector_indexing_suite.hpp>

#include "graph_filtering.hh"

#include "graph.hh"
#include "graph_selectors.hh"
#include "graph_properties.hh"
33
#include "graph_util.hh"
34
35
36

#include "random.hh"

37
#include "config.h"
38
39

#include "graph_blockmodel_overlap.hh"
40
41
42
#include "graph_blockmodel.hh"

using namespace boost;
43
using namespace graph_tool;
44

45
46


47
48
49
50
// ====================
// Entropy calculation
// ====================

51
52
53
54
55
56
57

// Repeated computation of x*log(x) and log(x) actually adds up to a lot of
// time. A significant speedup can be made by caching pre-computed values. This
// is doable since the values of mrse are bounded in [0, 2E], where E is the
// total number of edges in the network. Here we simply grow the cache as
// needed.

58
59
60
namespace graph_tool
{

61
62
63
64
65
66
67
68
69
70
71
72
73
vector<double> __safelog_cache;
vector<double> __xlogx_cache;
vector<double> __lgamma_cache;

void init_safelog(size_t x)
{
    size_t old_size = __safelog_cache.size();
    if (x >= old_size)
    {
        __safelog_cache.resize(x + 1);
        for (size_t i = old_size; i < __safelog_cache.size(); ++i)
            __safelog_cache[i] = safelog(double(i));
    }
74
75
}

76
void clear_safelog()
77
{
78
    vector<double>().swap(__safelog_cache);
79
80
}

81
82

void init_xlogx(size_t x)
83
{
84
85
86
87
88
89
90
    size_t old_size = __xlogx_cache.size();
    if (x >= old_size)
    {
        __xlogx_cache.resize(x + 1);
        for (size_t i = old_size; i < __xlogx_cache.size(); ++i)
            __xlogx_cache[i] = i * safelog(i);
    }
91
92
}

93
94
95
96
void clear_xlogx()
{
    vector<double>().swap(__xlogx_cache);
}
97

98
99
100
101
102
103
104
105
106
107
void init_lgamma(size_t x)
{
    size_t old_size = __lgamma_cache.size();
    if (x >= old_size)
    {
        __lgamma_cache.resize(x + 1);
        for (size_t i = old_size; i < __lgamma_cache.size(); ++i)
            __lgamma_cache[i] = lgamma(i);
    }
}
108

109
110
111
112
void clear_lgamma()
{
    vector<double>().swap(__lgamma_cache);
}
113

114
115
116
117
118

}

double do_get_ent(GraphInterface& gi, boost::any omrs, boost::any omrp,
                  boost::any omrm, boost::any owr, bool deg_corr)
119
120
121
122
123
124
125
126
127
128
129
130
{
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::vertex_index_map_t>::type
        vmap_t;
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::edge_index_map_t>::type
        emap_t;
    emap_t mrs = any_cast<emap_t>(omrs);
    vmap_t mrp = any_cast<vmap_t>(omrp);
    vmap_t mrm = any_cast<vmap_t>(omrm);
    vmap_t wr = any_cast<vmap_t>(owr);

131
132
    double S = 0;
    run_action<>()
Tiago Peixoto's avatar
Tiago Peixoto committed
133
134
135
        (gi, std::bind(entropy(), mrs, mrp, mrm, wr, deg_corr,
                       placeholders::_1,
                       std::ref(S)))();
136
    return S;
137
138
}

139
140
double do_get_ent_dense(GraphInterface& gi, boost::any omrs, boost::any owr,
                        bool multigraph)
141
142
143
144
145
146
147
148
149
150
{
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::vertex_index_map_t>::type
        vmap_t;
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::edge_index_map_t>::type
        emap_t;
    emap_t mrs = any_cast<emap_t>(omrs);
    vmap_t wr = any_cast<vmap_t>(owr);

151
152
    double S = 0;
    run_action<>()
Tiago Peixoto's avatar
Tiago Peixoto committed
153
154
        (gi, std::bind(entropy_dense(), mrs, wr, multigraph,
                       placeholders::_1, std::ref(S)))();
155
    return S;
156
157
}

158
159
160
161
162
163
164
165
166
167
168
169
170
double do_get_ent_parallel(GraphInterface& gi, boost::any oweight)
{
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::edge_index_map_t>::type
        emap_t;
    emap_t weight = any_cast<emap_t>(oweight);

    double S = 0;
    run_action<>()
        (gi, std::bind(entropy_parallel_edges(),
                       placeholders::_1, weight, std::ref(S)))();
    return S;
}
171

172

173
174
boost::any do_create_emat(GraphInterface& gi, boost::any ob,
                          GraphInterface& bgi)
175
176
{
    boost::any emat;
177
178
179
180
181
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::vertex_index_map_t>::type
        vmap_t;
    vmap_t b = any_cast<vmap_t>(ob);

182
    if (gi.get_directed())
183
184
185
    {
        run_action<>()(gi, std::bind<void>(create_emat(), placeholders::_1,
                                           std::ref(b),
186
                                           std::ref(bgi.get_graph()),
187
188
189
190
                                           std::ref(emat)))();
    }
    else
    {
191
        UndirectedAdaptor<GraphInterface::multigraph_t> ug(bgi.get_graph());
192
193
194
195
196
        run_action<>()(gi, std::bind<void>(create_emat(), placeholders::_1,
                                           std::ref(b),
                                           std::ref(ug),
                                           std::ref(emat)))();
    }
197
198
199
    return emat;
}

200
201
boost::any do_create_ehash(GraphInterface& gi, boost::any ob,
                           GraphInterface& bgi, rng_t& rng)
202
203
{
    boost::any emat;
204
205
206
207
208
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::vertex_index_map_t>::type
        vmap_t;
    vmap_t b = any_cast<vmap_t>(ob);

209
    if (gi.get_directed())
210
211
212
    {
        run_action<>()(gi, std::bind<void>(create_ehash(), placeholders::_1,
                                           std::ref(b),
213
                                           std::ref(bgi.get_graph()),
214
215
216
217
218
                                           std::ref(emat),
                                           std::ref(rng)))();
    }
    else
    {
219
        UndirectedAdaptor<GraphInterface::multigraph_t> ug(bgi.get_graph());
220
221
222
223
224
225
        run_action<>()(gi, std::bind<void>(create_ehash(), placeholders::_1,
                                           std::ref(b),
                                           std::ref(ug),
                                           std::ref(emat),
                                           std::ref(rng)))();
    }
226
227
    return emat;
}
228
229
230
231
232
233
234
235
236
237

//============
// Main loops
//============


template <class Eprop, class Vprop, class VEprop>
struct move_sweep_dispatch
{
    move_sweep_dispatch(Eprop eweight, Vprop vweight, boost::any egroups,
238
                        VEprop esrcpos, VEprop etgtpos, Vprop label,
239
240
241
242
                        vector<int64_t>& vlist, vector<int64_t>& block_list,
                        vector<int64_t>& target_list, bool deg_corr, bool dense,
                        bool multigraph, double beta, bool sequential,
                        bool parallel, bool random_move, double c, bool verbose,
243
                        size_t edge_index_range, size_t nmerges, size_t niter,
244
245
246
                        Vprop merge_map, partition_stats_t& partition_stats,
                        rng_t& rng, double& S, size_t& nmoves,
                        GraphInterface& bgi)
247
248

        : eweight(eweight), vweight(vweight), oegroups(egroups), esrcpos(esrcpos),
249
250
251
252
          etgtpos(etgtpos), label(label), vlist(vlist), block_list(block_list),
          target_list(target_list), deg_corr(deg_corr), dense(dense),
          multigraph(multigraph), beta(beta), sequential(sequential),
          parallel(parallel), random_move(random_move),
253
          c(c), verbose(verbose), edge_index_range(edge_index_range),
254
          nmerges(nmerges), niter(niter), merge_map(merge_map),
255
          partition_stats(partition_stats), rng(rng), S(S),
256
          nmoves(nmoves), bgi(bgi)
257
258
259
260
261
262
263
264
265
    {}

    Eprop eweight;
    Vprop vweight;
    boost::any oegroups;
    VEprop esrcpos;
    VEprop etgtpos;
    Vprop label;
    size_t n;
266
    vector<int64_t>& vlist;
267
    vector<int64_t>& block_list;
268
    vector<int64_t>& target_list;
269
    bool deg_corr;
270
271
    bool dense;
    bool multigraph;
272
273
    double beta;
    bool sequential;
274
    bool parallel;
275
    bool random_move;
276
    double c;
277
    bool verbose;
278
    size_t edge_index_range;
279
    size_t nmerges;
280
    size_t niter;
281
282
    Vprop merge_map;
    partition_stats_t& partition_stats;
283
284
285
286
287
288
289
    rng_t& rng;
    double& S;
    size_t& nmoves;
    GraphInterface& bgi;

    template <class Graph>
    void operator()(Eprop mrs, Vprop mrp, Vprop mrm, Vprop wr, Vprop b,
290
                    Graph& g, boost::any& emat, boost::any sampler,
291
                    boost::any cavity_sampler, bool weighted) const
292
293
294
    {
        if (is_directed::apply<Graph>::type::value)
        {
295
            dispatch(mrs, mrp, mrm, wr, b, g, emat, sampler, cavity_sampler,
296
                     bgi.get_graph(), weighted);
297
298
299
        }
        else
        {
300
            UndirectedAdaptor<GraphInterface::multigraph_t> ug(bgi.get_graph());
301
302
            dispatch(mrs, mrp, mrm, wr, b, g, emat, sampler, cavity_sampler, ug,
                     weighted);
303
304
305
306
307
        }
    }

    template <class Graph, class BGraph>
    void dispatch(Eprop mrs, Vprop mrp, Vprop mrm, Vprop wr, Vprop b, Graph& g,
308
309
                  boost::any& aemat, boost::any asampler,
                  boost::any acavity_sampler, BGraph& bg, bool weighted) const
310
311
312
313
314
315
    {
        if (weighted)
        {
            typedef typename property_map_type::apply<DynamicSampler<std::tuple<typename graph_traits<Graph>::edge_descriptor, bool> >,
                                                      GraphInterface::vertex_index_map_t>::type vemap_t;
            vemap_t egroups = any_cast<vemap_t>(oegroups);
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334

            try
            {
                typedef typename get_emat_t::apply<BGraph>::type emat_t;
                emat_t& emat = any_cast<emat_t&>(aemat);
                size_t B = num_vertices(bg);
                size_t max_BE = is_directed::apply<Graph>::type::value ?
                    B * B : (B * (B + 1)) / 2;
                dispatch(mrs.get_unchecked(max_BE), mrp, mrm, wr, b, g,
                         asampler, acavity_sampler, bg, egroups, emat);
            }
            catch (bad_any_cast&)
            {
                typedef typename get_ehash_t::apply<BGraph>::type emat_t;
                emat_t& emat = any_cast<emat_t&>(aemat);

                dispatch(mrs, mrp, mrm, wr, b, g, asampler, acavity_sampler, bg,
                         egroups, emat);
            }
335
336
337
338
339
340
        }
        else
        {
            typedef typename property_map_type::apply<vector<std::tuple<typename graph_traits<Graph>::edge_descriptor, bool> >,
                                                      GraphInterface::vertex_index_map_t>::type vemap_t;
            vemap_t egroups = any_cast<vemap_t>(oegroups);
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359

            try
            {
                typedef typename get_emat_t::apply<BGraph>::type emat_t;
                emat_t& emat = any_cast<emat_t&>(aemat);
                size_t B = num_vertices(bg);
                size_t max_BE = is_directed::apply<Graph>::type::value ?
                    B * B : (B * (B + 1)) / 2;
                dispatch(mrs.get_unchecked(max_BE), mrp, mrm, wr, b, g,
                         asampler, acavity_sampler, bg, egroups, emat);
            }
            catch (bad_any_cast&)
            {
                typedef typename get_ehash_t::apply<BGraph>::type emat_t;
                emat_t& emat = any_cast<emat_t&>(aemat);

                dispatch(mrs, mrp, mrm, wr, b, g, asampler, acavity_sampler, bg,
                         egroups, emat);
            }
360
361
362
        }
    }

363
364
365
366
    template <class Graph, class BGraph, class Egroups, class Emat, class MEprop>
    void dispatch(MEprop mrs, Vprop mrp, Vprop mrm, Vprop wr, Vprop b, Graph& g,
                  boost::any asampler, boost::any acavity_sampler, BGraph& bg,
                  Egroups egroups, Emat& emat) const
367
    {
368
369
        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;

370
        size_t eidx = random_move ? 1 : edge_index_range;
371
372

        typedef typename property_map<Graph, vertex_index_t>::type vindex_map_t;
Tiago Peixoto's avatar
Tiago Peixoto committed
373
        typedef typename property_map_type::apply<Sampler<vertex_t, boost::mpl::false_>,
374
375
376
                                                  vindex_map_t>::type::unchecked_t
            sampler_map_t;
        sampler_map_t sampler = any_cast<sampler_map_t>(asampler);
377
        sampler_map_t cavity_sampler = any_cast<sampler_map_t>(acavity_sampler);
378

379
380
        ConstantPropertyMap<std::array<int, 1>, typename graph_traits<Graph>::vertex_descriptor>
            cv(std::array<int, 1>({{-1}}));
381
382
383
384
385
386
387
388
        IdentityArrayPropertyMap<typename graph_traits<Graph>::vertex_descriptor> vmap;
        boost::typed_identity_property_map<int> identity;

        // make sure the properties are _unchecked_, since otherwise it
        // affects performance

        overlap_stats_t ostats;
        vector<size_t> free_blocks;
389
        auto state = make_block_state(g, eweight.get_unchecked(edge_index_range),
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
                                      vweight.get_unchecked(num_vertices(g)),
                                      b.get_unchecked(num_vertices(g)), bg,
                                      emat, mrs,
                                      mrp.get_unchecked(num_vertices(bg)),
                                      mrm.get_unchecked(num_vertices(bg)),
                                      wr.get_unchecked(num_vertices(bg)),
                                      egroups.get_unchecked(num_vertices(bg)),
                                      esrcpos.get_unchecked(eidx),
                                      etgtpos.get_unchecked(eidx), sampler,
                                      cavity_sampler, partition_stats, ostats,
                                      identity, identity, free_blocks,
                                      false, false, true);

        vector<decltype(state)> states = {state};
        vector<EntrySet<Graph>> m_entries = {EntrySet<Graph>(num_vertices(bg))};

406
407
408
409
410
411
        if (target_list.empty())
        {
            move_sweep(states, m_entries,
                       wr.get_unchecked(num_vertices(bg)),
                       b.get_unchecked(num_vertices(g)),
                       cv, vmap,
412
413
414
                       label.get_unchecked(num_vertices(bg)),
                       vlist, block_list,
                       deg_corr, dense, multigraph, beta,
415
                       eweight.get_unchecked(edge_index_range),
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
                       vweight.get_unchecked(num_vertices(g)),
                       g, sequential, parallel, random_move, c,
                       nmerges,
                       merge_map.get_unchecked(num_vertices(g)),
                       niter, num_vertices(bg),
                       verbose, rng, S, nmoves);
        }
        else
        {
            auto ub = b.get_unchecked(num_vertices(g));
            for (size_t i = 0; i < vlist.size(); ++i)
            {
                size_t v = vlist[i];
                size_t s = target_list[i];
                S += virtual_move(v, s, ub, cv, vmap, states, m_entries,
                                  dense, deg_corr, multigraph);
                move_vertex(v, s, ub, cv, vmap, deg_corr, states,
                            not random_move);
                nmoves++;
            }
        }
437
438
439
440
    }
};


Tiago Peixoto's avatar
Tiago Peixoto committed
441
442
boost::python::object do_move_sweep(GraphInterface& gi, GraphInterface& bgi,
                                    boost::any& emat, boost::any sampler,
443
444
445
                                    boost::any cavity_sampler, boost::any omrs,
                                    boost::any omrp, boost::any omrm,
                                    boost::any owr, boost::any ob,
446
447
448
                                    boost::any olabel,
                                    vector<int64_t>& vlist,
                                    vector<int64_t>& block_list,
449
                                    vector<int64_t>& target_list,
450
451
452
453
454
455
                                    bool deg_corr, bool dense, bool multigraph,
                                    boost::any oeweight, boost::any ovweight,
                                    boost::any oegroups, boost::any oesrcpos,
                                    boost::any oetgtpos, double beta,
                                    bool sequential, bool parallel,
                                    bool random_move, double c, bool weighted,
456
457
                                    size_t nmerges, boost::any omerge_map,
                                    size_t niter,
458
459
                                    partition_stats_t& partition_stats,
                                    bool verbose, rng_t& rng)
460
461
462
463
464
465
466
{
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::vertex_index_map_t>::type
        vmap_t;
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::edge_index_map_t>::type
        emap_t;
467
    typedef property_map_type::apply<int32_t,
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
                                     GraphInterface::edge_index_map_t>::type
        vemap_t;
    emap_t mrs = any_cast<emap_t>(omrs);
    vmap_t mrp = any_cast<vmap_t>(omrp);
    vmap_t mrm = any_cast<vmap_t>(omrm);
    vmap_t wr = any_cast<vmap_t>(owr);
    vmap_t b = any_cast<vmap_t>(ob);
    vmap_t label = any_cast<vmap_t>(olabel);
    emap_t eweight = any_cast<emap_t>(oeweight);
    vmap_t vweight = any_cast<vmap_t>(ovweight);

    vemap_t esrcpos = any_cast<vemap_t>(oesrcpos);
    vemap_t etgtpos = any_cast<vemap_t>(oetgtpos);

    double S = 0;
    size_t nmoves = 0;
484

485
486
    vmap_t merge_map = any_cast<vmap_t>(omerge_map);

Tiago Peixoto's avatar
Tiago Peixoto committed
487
488
489
    run_action<graph_tool::detail::all_graph_views, boost::mpl::true_>()
        (gi, std::bind(move_sweep_dispatch<emap_t, vmap_t, vemap_t>
                       (eweight, vweight, oegroups, esrcpos, etgtpos,
490
491
                        label, vlist, block_list,
                        target_list, deg_corr, dense, multigraph,
492
                        beta, sequential, parallel, random_move, c, verbose,
493
                        gi.get_edge_index_range(), nmerges, niter, merge_map,
494
                        partition_stats, rng, S, nmoves, bgi),
Tiago Peixoto's avatar
Tiago Peixoto committed
495
                       mrs, mrp, mrm, wr, b, placeholders::_1,
496
                       std::ref(emat), sampler, cavity_sampler, weighted))();
Tiago Peixoto's avatar
Tiago Peixoto committed
497
    return boost::python::make_tuple(S, nmoves);
498
499
}

500
501
502
503
504
struct build_egroups
{
    template <class Eprop, class Vprop, class VEprop, class Graph, class VertexIndex>
    void operator()(Vprop b, boost::any& oegroups, VEprop esrcpos,
                    VEprop etgtpos, Eprop eweight, Graph& g,
505
                    VertexIndex vertex_index, size_t B, bool weighted, bool empty) const
506
507
    {
        egroups_manage::build(b, oegroups, esrcpos, etgtpos, eweight, g,
508
                              vertex_index, B, weighted, empty);
509
510
511
    }
};

512
513
boost::any do_build_egroups(GraphInterface& gi, GraphInterface& bgi,
                            boost::any ob, boost::any oeweights,
514
                            boost::any oesrcpos, boost::any oetgtpos,
515
                            bool weighted, bool empty)
516
517
518
519
520
521
522
{
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::vertex_index_map_t>::type
        vmap_t;
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::edge_index_map_t>::type
        emap_t;
523
    typedef property_map_type::apply<int32_t,
524
525
526
527
528
529
530
531
532
                                     GraphInterface::edge_index_map_t>::type
        vemap_t;
    vmap_t b = any_cast<vmap_t>(ob);

    vemap_t esrcpos = any_cast<vemap_t>(oesrcpos);
    vemap_t etgtpos = any_cast<vemap_t>(oetgtpos);
    emap_t eweights = any_cast<emap_t>(oeweights);

    boost::any oegroups;
Tiago Peixoto's avatar
Tiago Peixoto committed
533
534
    run_action<graph_tool::detail::all_graph_views, boost::mpl::true_>()
        (gi, std::bind<void>(build_egroups(), b, std::ref(oegroups),
535
536
537
                             esrcpos.get_unchecked(gi.get_edge_index_range()),
                             etgtpos.get_unchecked(gi.get_edge_index_range()),
                             eweights.get_unchecked(gi.get_edge_index_range()),
538
539
                             placeholders::_1, bgi.get_vertex_index(),
                             bgi.get_num_vertices(), weighted, empty))();
540
541
542
    return oegroups;
}

543
boost::any do_init_neighbour_sampler(GraphInterface& gi, boost::any oeweights,
544
                                     bool self_loops, bool empty)
545
546
547
548
549
550
551
{
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::edge_index_map_t>::type
        emap_t;
    emap_t eweights = any_cast<emap_t>(oeweights);

    boost::any osampler;
Tiago Peixoto's avatar
Tiago Peixoto committed
552
553
    run_action<graph_tool::detail::all_graph_views, boost::mpl::true_>()
        (gi, std::bind(init_neighbour_sampler(), placeholders::_1, eweights,
554
                       self_loops, empty, std::ref(osampler)))();
555
556
557
    return osampler;
}

558
559
560
561
struct collect_edge_marginals_dispatch
{
    template <class Graph, class Vprop, class MEprop>
    void operator()(Graph& g, size_t B, Vprop cb, MEprop p,
Tiago Peixoto's avatar
Tiago Peixoto committed
562
                    std::tuple<boost::any, GraphInterface&> abg) const
563
    {
564
        Graph& bg = *any_cast<Graph*>(get<0>(abg));
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
        collect_edge_marginals(B, cb.get_unchecked(num_vertices(bg)), p, g, bg);
    }
};

void do_collect_edge_marginals(GraphInterface& gi, GraphInterface& gbi,
                               size_t B, boost::any ob, boost::any op)
{
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::vertex_index_map_t>::type
        vmap_t;
    typedef property_map_type::apply<vector<int32_t>,
                                     GraphInterface::edge_index_map_t>::type
        emap_t;
    vmap_t b = any_cast<vmap_t>(ob);
    emap_t p = any_cast<emap_t>(op);

Tiago Peixoto's avatar
Tiago Peixoto committed
581
582
583
    run_action<graph_tool::detail::all_graph_views, boost::mpl::true_>()
        (gi, std::bind<void>(collect_edge_marginals_dispatch(),
                             placeholders::_1, B, b, p,
584
                             std::tuple<boost::any, GraphInterface&>(gbi.get_graph_view(), gbi)))();
585
586
}

Tiago Peixoto's avatar
Tiago Peixoto committed
587
588
boost::python::tuple do_bethe_entropy(GraphInterface& gi, size_t B, boost::any op,
                                      boost::any opv)
589
590
591
592
593
594
595
596
597
598
599
{
    typedef property_map_type::apply<vector<double>,
                                     GraphInterface::vertex_index_map_t>::type
        vmap_t;
    typedef property_map_type::apply<vector<int32_t>,
                                     GraphInterface::edge_index_map_t>::type
        emap_t;
    emap_t p = any_cast<emap_t>(op);
    vmap_t pv = any_cast<vmap_t>(opv);

    double H=0, sH=0, Hmf=0, sHmf=0;
Tiago Peixoto's avatar
Tiago Peixoto committed
600
601
602
603
604
    run_action<graph_tool::detail::all_graph_views, boost::mpl::true_>()
        (gi, std::bind<void>(bethe_entropy(),
                             placeholders::_1, B, p, pv, std::ref(H), std::ref(sH),
                             std::ref(Hmf), std::ref(sHmf)))();
    return boost::python::make_tuple(H, sH, Hmf, sHmf);
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
}


struct collect_vertex_marginals_dispatch
{
    template <class Graph, class Vprop, class MEprop>
    void operator()(Graph& g, Vprop cb, MEprop p) const
    {
        collect_vertex_marginals(cb.get_unchecked(num_vertices(g)),
                                 p.get_unchecked(num_vertices(g)), g);
    }
};


void do_collect_vertex_marginals(GraphInterface& gi, boost::any ob,
                                 boost::any op)
{
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::vertex_index_map_t>::type
        vmap_t;
    vmap_t b = any_cast<vmap_t>(ob);

Tiago Peixoto's avatar
Tiago Peixoto committed
627
628
629
    run_action<graph_tool::detail::all_graph_views, boost::mpl::true_>()
        (gi, std::bind(collect_vertex_marginals_dispatch(),
                       placeholders::_1, b, placeholders::_2),
630
631
632
         vertex_scalar_vector_properties())(op);
}

633
634

struct get_deg_entropy_term
635
{
636
637
    template <class Graph, class Emap, class VImap>
    void operator()(Graph& g, Emap eweight, VImap ignore_degrees, double& S) const
638
639
640
    {
        for(auto v : vertices_range(g))
        {
641
            if (ignore_degrees[v] == 1)
642
643
                continue;
            S -= lgamma_fast(in_degreeS()(v, g, eweight) + 1);
644
645
646
            if (ignore_degrees[v] == 2)
                continue;
            S -= lgamma_fast(out_degreeS()(v, g, eweight) + 1);
647
648
649
650
651
652
653
        }
    }
};

struct get_deg_entropy_term_overlap
{

654
    template <class Graph, class Vprop>
655
656
    void operator()(Graph& g, Vprop b, overlap_stats_t& overlap_stats, size_t N,
                    double& S) const
657
    {
658
        typedef gt_hash_map<int, int> map_t;
659
660

        map_t in_hist, out_hist;
661

662
        for (size_t v = 0; v < N; ++v)
663
        {
664
665
            in_hist.clear();
            out_hist.clear();
666

667
668
            auto& half_edges = overlap_stats.get_half_edges(v);
            for (size_t u : half_edges)
669
            {
670
671
                in_hist[b[u]] += in_degreeS()(u, g);
                out_hist[b[u]] += out_degree(u, g);
672
            }
673
674
675
676
677

            for (auto& k_c : in_hist)
                S -= lgamma_fast(k_c.second + 1);
            for (auto& k_c : out_hist)
                S -= lgamma_fast(k_c.second + 1);
678
679
680
681
        }
    }
};

682
double do_get_deg_entropy_term(GraphInterface& gi, boost::any ob,
683
684
                               overlap_stats_t& overlap_stats, size_t N,
                               boost::any aeweight, boost::any aignore_degrees)
685
686
687
688
{
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::vertex_index_map_t>::type
        vmap_t;
689
690
691
692
693
694
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::edge_index_map_t>::type
        emap_t;
    typedef property_map_type::apply<uint8_t,
                                     GraphInterface::vertex_index_map_t>::type
        vimap_t;
695
696

    double S = 0;
697
698
699
700
701
702
703
704
705
706
707

    if (!ob.empty())
    {
        vmap_t b = any_cast<vmap_t>(ob);
        run_action<>()
            (gi, std::bind(get_deg_entropy_term_overlap(),
                           placeholders::_1, b, std::ref(overlap_stats), N,
                           std::ref(S)))();
    }
    else
    {
708
709
        emap_t eweight = any_cast<emap_t>(aeweight);
        vimap_t ignore_degrees = any_cast<vimap_t>(aignore_degrees);
710
711
        run_action<>()
            (gi, std::bind(get_deg_entropy_term(),
712
713
                           placeholders::_1, eweight, ignore_degrees,
                           std::ref(S)))();
714
715
    }

716
717
718
719
    return S;
}


720
vector<int64_t> get_vector(size_t n)
721
{
722
    return vector<int64_t>(n);
723
724
}

725
template <class Value>
Tiago Peixoto's avatar
Tiago Peixoto committed
726
void vector_map(boost::python::object ovals, boost::python::object omap)
727
728
{

729
730
    multi_array_ref<Value,1> vals = get_array<Value,1>(ovals);
    multi_array_ref<Value,1> map = get_array<Value,1>(omap);
731
732
733
734

    size_t pos = 0;
    for (size_t i = 0; i < vals.size(); ++i)
    {
735
        Value v = vals[i];
736
737
738
739
740
741
        if (map[v] == -1)
            map[v] = pos++;
        vals[i] = map[v];
    }
}

742
template <class Value>
743
744
745
void vector_continuous_map(boost::python::object ovals)
{

746
    multi_array_ref<Value,1> vals = get_array<Value,1>(ovals);
747
    gt_hash_map<Value, size_t> map;
748
749
750

    for (size_t i = 0; i < vals.size(); ++i)
    {
751
        Value v = vals[i];
752
753
754
755
756
757
758
        auto iter = map.find(v);
        if (iter == map.end())
            iter = map.insert(make_pair(v, map.size())).first;
        vals[i] = iter->second;
    }
}

759
template <class Value>
Tiago Peixoto's avatar
Tiago Peixoto committed
760
void vector_rmap(boost::python::object ovals, boost::python::object omap)
761
762
{

763
764
    multi_array_ref<Value,1> vals = get_array<Value,1>(ovals);
    multi_array_ref<Value,1> map = get_array<Value,1>(omap);
765
766
767
768
769
770
771

    for (size_t i = 0; i < vals.size(); ++i)
    {
        map[vals[i]] = i;
    }
}

772
773
774
775
776
777
778
779
780
python::tuple python_get_mu_l(double N, double E, double epsilon=1e-8)
{
    double mu, l;
    get_mu_l(N, E, mu, l, epsilon);
    return python::make_tuple(mu, l);
}

struct get_partition_stats
{
781
    template <class Graph, class Vprop, class Eprop, class Mprop>
782
    void operator()(Graph& g, Vprop b, Eprop eweight, size_t N, size_t B,
783
784
                    bool edges_dl, partition_stats_t& partition_stats,
                    Mprop ignore_degrees) const
785
    {
786
787
        partition_stats = partition_stats_t(g, b, eweight, N, B, edges_dl,
                                            ignore_degrees);
788
789
790
791
792
    }
};

partition_stats_t
do_get_partition_stats(GraphInterface& gi, boost::any ob, boost::any aeweight,
793
794
                       size_t N, size_t B, bool edges_dl,
                       boost::any aignore_degrees)
795
796
797
798
799
800
801
{
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::vertex_index_map_t>::type
        vmap_t;
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::edge_index_map_t>::type
        emap_t;
802
803
804
    typedef property_map_type::apply<uint8_t,
                                     GraphInterface::vertex_index_map_t>::type
        mvmap_t;
805
806
807
808
809

    partition_stats_t partition_stats;

    vmap_t b = any_cast<vmap_t>(ob);
    emap_t eweight = any_cast<emap_t>(aeweight);
810
    mvmap_t ignore_degrees = any_cast<mvmap_t>(aignore_degrees);
811
812

    run_action<>()(gi, std::bind(get_partition_stats(),
813
                                 placeholders::_1, b, eweight, N, B, edges_dl,
814
                                 std::ref(partition_stats), ignore_degrees))();
815
816
817
    return partition_stats;
}

818
819
820
821
void export_blockmodel()
{
    using namespace boost::python;

822
823
824
825
826
827
828
    class_<partition_stats_t>("partition_stats")
        .def("is_enabled", &partition_stats_t::is_enabled)
        .def("get_partition_dl", &partition_stats_t::get_partition_dl)
        .def("get_deg_dl", &partition_stats_t::get_deg_dl);

    def("init_partition_stats", do_get_partition_stats);

829
830
831
832
833
834
    def("init_safelog", init_safelog);
    def("clear_safelog", clear_safelog);
    def("init_xlogx", init_xlogx);
    def("clear_xlogx", clear_xlogx);
    def("init_lgamma", init_lgamma);
    def("clear_lgamma", clear_lgamma);
835
836
837
838
    def("get_xi", get_xi<double,double>);
    def("get_xi_fast", get_xi_fast<double,double>);
    def("get_mu_l", python_get_mu_l);
    def("polylog", polylog<double>);
839
840
841
    def("lbinom_careful", lbinom_careful);
    def("lbinom_fast", lbinom_fast);
    def("lbinom", lbinom);
842

843
    def("get_vector", get_vector);
844
845
846
847
848
849
    def("vector_map", vector_map<int32_t>);
    def("vector_map64", vector_map<int64_t>);
    def("vector_rmap", vector_rmap<int32_t>);
    def("vector_rmap64", vector_rmap<int64_t>);
    def("vector_continuous_map", vector_continuous_map<int32_t>);
    def("vector_continuous_map64", vector_continuous_map<int64_t>);
850
851

    def("create_emat", do_create_emat);
852
    def("create_ehash", do_create_ehash);
853
    def("build_egroups", do_build_egroups);
854
    def("init_neighbour_sampler", do_init_neighbour_sampler);
855
856
857
858

    def("move_sweep", do_move_sweep);

    def("entropy", do_get_ent);
859
    def("entropy_dense", do_get_ent_dense);
860
861
    def("entropy_parallel", do_get_ent_parallel);
    def("deg_entropy_term", do_get_deg_entropy_term);
862
863
864
865
866

    def("edge_marginals", do_collect_edge_marginals);
    def("bethe_entropy", do_bethe_entropy);
    def("vertex_marginals", do_collect_vertex_marginals);
}