graph_blockmodel.cc 28.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
boost::any do_create_emat(GraphInterface& gi)
174
175
{
    boost::any emat;
Tiago Peixoto's avatar
Tiago Peixoto committed
176
177
    run_action<>()(gi, std::bind<void>(create_emat(), placeholders::_1,
                                       std::ref(emat)))();
178
179
180
    return emat;
}

181
182
183
boost::any do_create_ehash(GraphInterface& gi)
{
    boost::any emat;
Tiago Peixoto's avatar
Tiago Peixoto committed
184
185
    run_action<>()(gi, std::bind<void>(create_ehash(), placeholders::_1,
                                       std::ref(emat)))();
186
187
    return emat;
}
188
189
190
191
192
193
194
195
196
197

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


template <class Eprop, class Vprop, class VEprop>
struct move_sweep_dispatch
{
    move_sweep_dispatch(Eprop eweight, Vprop vweight, boost::any egroups,
198
                        VEprop esrcpos, VEprop etgtpos, Vprop label,
199
200
201
202
203
204
205
                        vector<int64_t>& vlist, 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, size_t max_edge_index,
                        size_t nmerges, size_t niter, Vprop merge_map,
                        partition_stats_t& partition_stats, rng_t& rng,
                        double& S, size_t& nmoves, GraphInterface& bgi)
206
207

        : eweight(eweight), vweight(vweight), oegroups(egroups), esrcpos(esrcpos),
208
          etgtpos(etgtpos), label(label), vlist(vlist), target_list(target_list), 
209
          deg_corr(deg_corr), dense(dense), multigraph(multigraph), beta(beta),
210
211
          sequential(sequential), parallel(parallel), random_move(random_move),
          c(c), verbose(verbose), max_edge_index(max_edge_index),
212
          nmerges(nmerges), niter(niter), merge_map(merge_map),
213
          partition_stats(partition_stats), rng(rng), S(S),
214
          nmoves(nmoves), bgi(bgi)
215
216
217
218
219
220
221
222
223
    {}

    Eprop eweight;
    Vprop vweight;
    boost::any oegroups;
    VEprop esrcpos;
    VEprop etgtpos;
    Vprop label;
    size_t n;
224
225
    vector<int64_t>& vlist;
    vector<int64_t>& target_list;
226
    bool deg_corr;
227
228
    bool dense;
    bool multigraph;
229
230
    double beta;
    bool sequential;
231
    bool parallel;
232
    bool random_move;
233
    double c;
234
235
    bool verbose;
    size_t max_edge_index;
236
    size_t nmerges;
237
    size_t niter;
238
239
    Vprop merge_map;
    partition_stats_t& partition_stats;
240
241
242
243
244
245
246
    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,
247
                    Graph& g, boost::any& emat, boost::any sampler,
248
                    boost::any cavity_sampler, bool weighted) const
249
250
251
    {
        if (is_directed::apply<Graph>::type::value)
        {
252
253
            dispatch(mrs, mrp, mrm, wr, b, g, emat, sampler, cavity_sampler,
                     bgi.GetGraph(), weighted);
254
255
256
257
        }
        else
        {
            UndirectedAdaptor<GraphInterface::multigraph_t> ug(bgi.GetGraph());
258
259
            dispatch(mrs, mrp, mrm, wr, b, g, emat, sampler, cavity_sampler, ug,
                     weighted);
260
261
262
263
264
        }
    }

    template <class Graph, class BGraph>
    void dispatch(Eprop mrs, Vprop mrp, Vprop mrm, Vprop wr, Vprop b, Graph& g,
265
266
                  boost::any& aemat, boost::any asampler,
                  boost::any acavity_sampler, BGraph& bg, bool weighted) const
267
268
269
270
271
272
    {
        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);
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291

            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);
            }
292
293
294
295
296
297
        }
        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);
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316

            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);
            }
317
318
319
        }
    }

320
321
322
323
    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
324
    {
325
326
327
328
329
        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;

        size_t eidx = random_move ? 1 : max_edge_index;

        typedef typename property_map<Graph, vertex_index_t>::type vindex_map_t;
Tiago Peixoto's avatar
Tiago Peixoto committed
330
        typedef typename property_map_type::apply<Sampler<vertex_t, boost::mpl::false_>,
331
332
333
                                                  vindex_map_t>::type::unchecked_t
            sampler_map_t;
        sampler_map_t sampler = any_cast<sampler_map_t>(asampler);
334
        sampler_map_t cavity_sampler = any_cast<sampler_map_t>(acavity_sampler);
335

336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
        ConstantPropertyMap<std::array<int, 1>, typename graph_traits<Graph>::vertex_descriptor> cv({-1});
        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;
        auto state = make_block_state(g, eweight.get_unchecked(max_edge_index),
                                      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))};

362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
        if (target_list.empty())
        {
            move_sweep(states, m_entries,
                       wr.get_unchecked(num_vertices(bg)),
                       b.get_unchecked(num_vertices(g)),
                       cv, vmap,
                       label.get_unchecked(num_vertices(bg)), vlist, deg_corr,
                       dense, multigraph, beta,
                       eweight.get_unchecked(max_edge_index),
                       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++;
            }
        }
392
393
394
395
    }
};


Tiago Peixoto's avatar
Tiago Peixoto committed
396
397
boost::python::object do_move_sweep(GraphInterface& gi, GraphInterface& bgi,
                                    boost::any& emat, boost::any sampler,
398
399
400
                                    boost::any cavity_sampler, boost::any omrs,
                                    boost::any omrp, boost::any omrm,
                                    boost::any owr, boost::any ob,
401
402
                                    boost::any olabel, vector<int64_t>& vlist,
                                    vector<int64_t>& target_list,
403
404
405
406
407
408
                                    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,
409
410
                                    size_t nmerges, boost::any omerge_map,
                                    size_t niter,
411
412
                                    partition_stats_t& partition_stats,
                                    bool verbose, rng_t& rng)
413
414
415
416
417
418
419
{
    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;
420
    typedef property_map_type::apply<int32_t,
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
                                     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;
437

438
439
    vmap_t merge_map = any_cast<vmap_t>(omerge_map);

Tiago Peixoto's avatar
Tiago Peixoto committed
440
441
442
    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,
443
444
                        label, vlist,  target_list, deg_corr, dense, multigraph,
                        beta, sequential, parallel, random_move, c, verbose,
445
                        gi.GetMaxEdgeIndex(), nmerges, niter, merge_map,
446
                        partition_stats, rng, S, nmoves, bgi),
Tiago Peixoto's avatar
Tiago Peixoto committed
447
                       mrs, mrp, mrm, wr, b, placeholders::_1,
448
                       std::ref(emat), sampler, cavity_sampler, weighted))();
Tiago Peixoto's avatar
Tiago Peixoto committed
449
    return boost::python::make_tuple(S, nmoves);
450
451
}

452
453
454
455
456
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,
457
                    VertexIndex vertex_index, size_t B, bool weighted, bool empty) const
458
459
    {
        egroups_manage::build(b, oegroups, esrcpos, etgtpos, eweight, g,
460
                              vertex_index, B, weighted, empty);
461
462
463
    }
};

464
465
boost::any do_build_egroups(GraphInterface& gi, GraphInterface& bgi,
                            boost::any ob, boost::any oeweights,
466
                            boost::any oesrcpos, boost::any oetgtpos,
467
                            bool weighted, bool empty)
468
469
470
471
472
473
474
{
    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;
475
    typedef property_map_type::apply<int32_t,
476
477
478
479
480
481
482
483
484
                                     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
485
486
487
488
489
    run_action<graph_tool::detail::all_graph_views, boost::mpl::true_>()
        (gi, std::bind<void>(build_egroups(), b, std::ref(oegroups),
                             esrcpos.get_unchecked(gi.GetMaxEdgeIndex()),
                             etgtpos.get_unchecked(gi.GetMaxEdgeIndex()),
                             eweights.get_unchecked(gi.GetMaxEdgeIndex()),
490
491
                             placeholders::_1, bgi.GetVertexIndex(),
                             bgi.GetNumberOfVertices(), weighted, empty))();
492
493
494
    return oegroups;
}

495
boost::any do_init_neighbour_sampler(GraphInterface& gi, boost::any oeweights,
496
                                     bool self_loops, bool empty)
497
498
499
500
501
502
503
{
    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
504
505
    run_action<graph_tool::detail::all_graph_views, boost::mpl::true_>()
        (gi, std::bind(init_neighbour_sampler(), placeholders::_1, eweights,
506
                       self_loops, empty, std::ref(osampler)))();
507
508
509
    return osampler;
}

510
511
512
513
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
514
                    std::tuple<boost::any, GraphInterface&> abg) const
515
    {
516
        Graph& bg = *any_cast<Graph*>(get<0>(abg));
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
        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
533
534
535
536
    run_action<graph_tool::detail::all_graph_views, boost::mpl::true_>()
        (gi, std::bind<void>(collect_edge_marginals_dispatch(),
                             placeholders::_1, B, b, p,
                             std::tuple<boost::any, GraphInterface&>(gbi.GetGraphView(), gbi)))();
537
538
}

Tiago Peixoto's avatar
Tiago Peixoto committed
539
540
boost::python::tuple do_bethe_entropy(GraphInterface& gi, size_t B, boost::any op,
                                      boost::any opv)
541
542
543
544
545
546
547
548
549
550
551
{
    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
552
553
554
555
556
    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);
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
}


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
579
580
581
    run_action<graph_tool::detail::all_graph_views, boost::mpl::true_>()
        (gi, std::bind(collect_vertex_marginals_dispatch(),
                       placeholders::_1, b, placeholders::_2),
582
583
584
         vertex_scalar_vector_properties())(op);
}

585
586

struct get_deg_entropy_term
587
{
588
589
    template <class Graph, class Emap, class VImap>
    void operator()(Graph& g, Emap eweight, VImap ignore_degrees, double& S) const
590
591
592
    {
        for(auto v : vertices_range(g))
        {
593
            if (ignore_degrees[v] == 1)
594
595
                continue;
            S -= lgamma_fast(in_degreeS()(v, g, eweight) + 1);
596
597
598
            if (ignore_degrees[v] == 2)
                continue;
            S -= lgamma_fast(out_degreeS()(v, g, eweight) + 1);
599
600
601
602
603
604
605
        }
    }
};

struct get_deg_entropy_term_overlap
{

606
    template <class Graph, class Vprop>
607
608
    void operator()(Graph& g, Vprop b, overlap_stats_t& overlap_stats, size_t N,
                    double& S) const
609
    {
610
        typedef gt_hash_map<int, int> map_t;
611
612

        map_t in_hist, out_hist;
613

614
        for (size_t v = 0; v < N; ++v)
615
        {
616
617
            in_hist.clear();
            out_hist.clear();
618

619
620
            auto& half_edges = overlap_stats.get_half_edges(v);
            for (size_t u : half_edges)
621
            {
622
623
                in_hist[b[u]] += in_degreeS()(u, g);
                out_hist[b[u]] += out_degree(u, g);
624
            }
625
626
627
628
629

            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);
630
631
632
633
        }
    }
};

634
double do_get_deg_entropy_term(GraphInterface& gi, boost::any ob,
635
636
                               overlap_stats_t& overlap_stats, size_t N,
                               boost::any aeweight, boost::any aignore_degrees)
637
638
639
640
{
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::vertex_index_map_t>::type
        vmap_t;
641
642
643
644
645
646
    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;
647
648

    double S = 0;
649
650
651
652
653
654
655
656
657
658
659

    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
    {
660
661
        emap_t eweight = any_cast<emap_t>(aeweight);
        vimap_t ignore_degrees = any_cast<vimap_t>(aignore_degrees);
662
663
        run_action<>()
            (gi, std::bind(get_deg_entropy_term(),
664
665
                           placeholders::_1, eweight, ignore_degrees,
                           std::ref(S)))();
666
667
    }

668
669
670
671
    return S;
}


672
vector<int64_t> get_vector(size_t n)
673
{
674
    return vector<int64_t>(n);
675
676
}

677
template <class Value>
Tiago Peixoto's avatar
Tiago Peixoto committed
678
void vector_map(boost::python::object ovals, boost::python::object omap)
679
680
{

681
682
    multi_array_ref<Value,1> vals = get_array<Value,1>(ovals);
    multi_array_ref<Value,1> map = get_array<Value,1>(omap);
683
684
685
686

    size_t pos = 0;
    for (size_t i = 0; i < vals.size(); ++i)
    {
687
        Value v = vals[i];
688
689
690
691
692
693
        if (map[v] == -1)
            map[v] = pos++;
        vals[i] = map[v];
    }
}

694
template <class Value>
695
696
697
void vector_continuous_map(boost::python::object ovals)
{

698
    multi_array_ref<Value,1> vals = get_array<Value,1>(ovals);
699
    gt_hash_map<Value, size_t> map;
700
701
702

    for (size_t i = 0; i < vals.size(); ++i)
    {
703
        Value v = vals[i];
704
705
706
707
708
709
710
        auto iter = map.find(v);
        if (iter == map.end())
            iter = map.insert(make_pair(v, map.size())).first;
        vals[i] = iter->second;
    }
}

711
template <class Value>
Tiago Peixoto's avatar
Tiago Peixoto committed
712
void vector_rmap(boost::python::object ovals, boost::python::object omap)
713
714
{

715
716
    multi_array_ref<Value,1> vals = get_array<Value,1>(ovals);
    multi_array_ref<Value,1> map = get_array<Value,1>(omap);
717
718
719
720
721
722
723

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

724
725
726
727
728
729
730
731
732
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
{
733
    template <class Graph, class Vprop, class Eprop, class Mprop>
734
    void operator()(Graph& g, Vprop b, Eprop eweight, size_t N, size_t B,
735
736
                    bool edges_dl, partition_stats_t& partition_stats,
                    Mprop ignore_degrees) const
737
    {
738
739
        partition_stats = partition_stats_t(g, b, eweight, N, B, edges_dl,
                                            ignore_degrees);
740
741
742
743
744
    }
};

partition_stats_t
do_get_partition_stats(GraphInterface& gi, boost::any ob, boost::any aeweight,
745
746
                       size_t N, size_t B, bool edges_dl,
                       boost::any aignore_degrees)
747
748
749
750
751
752
753
{
    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;
754
755
756
    typedef property_map_type::apply<uint8_t,
                                     GraphInterface::vertex_index_map_t>::type
        mvmap_t;
757
758
759
760
761

    partition_stats_t partition_stats;

    vmap_t b = any_cast<vmap_t>(ob);
    emap_t eweight = any_cast<emap_t>(aeweight);
762
    mvmap_t ignore_degrees = any_cast<mvmap_t>(aignore_degrees);
763
764

    run_action<>()(gi, std::bind(get_partition_stats(),
765
                                 placeholders::_1, b, eweight, N, B, edges_dl,
766
                                 std::ref(partition_stats), ignore_degrees))();
767
768
769
    return partition_stats;
}

770
771
772
773
void export_blockmodel()
{
    using namespace boost::python;

774
775
776
777
778
779
780
    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);

781
782
783
784
785
786
    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);
787
788
789
790
    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>);
791
792
793
    def("lbinom_careful", lbinom_careful);
    def("lbinom_fast", lbinom_fast);
    def("lbinom", lbinom);
794

795
    def("get_vector", get_vector);
796
797
798
799
800
801
    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>);
802
803

    def("create_emat", do_create_emat);
804
    def("create_ehash", do_create_ehash);
805
    def("build_egroups", do_build_egroups);
806
    def("init_neighbour_sampler", do_init_neighbour_sampler);
807
808
809
810

    def("move_sweep", do_move_sweep);

    def("entropy", do_get_ent);
811
    def("entropy_dense", do_get_ent_dense);
812
813
    def("entropy_parallel", do_get_ent_parallel);
    def("deg_entropy_term", do_get_deg_entropy_term);
814
815
816
817
818

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