graph_blockmodel.cc 27.1 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 199 200
                        VEprop esrcpos, VEprop etgtpos, Vprop label,
                        vector<int>& vlist, bool deg_corr, bool dense,
                        bool multigraph, double beta, bool sequential,
201
                        bool parallel, bool random_move, double c, bool verbose,
202
                        size_t max_edge_index, size_t nmerges, size_t niter,
203 204 205
                        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 209
          etgtpos(etgtpos), label(label), vlist(vlist),
          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 224 225
    {}

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

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

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

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

319 320 321 322
    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
323
    {
324 325 326 327 328
        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
329
        typedef typename property_map_type::apply<Sampler<vertex_t, boost::mpl::false_>,
330 331 332
                                                  vindex_map_t>::type::unchecked_t
            sampler_map_t;
        sampler_map_t sampler = any_cast<sampler_map_t>(asampler);
333
        sampler_map_t cavity_sampler = any_cast<sampler_map_t>(acavity_sampler);
334

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 362 363
        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))};

        move_sweep(states, m_entries,
                   wr.get_unchecked(num_vertices(bg)),
                   b.get_unchecked(num_vertices(g)),
364
                   cv, vmap,
365 366 367 368 369 370 371 372
                   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),
373
                   verbose, rng, S, nmoves);
374 375 376 377
    }
};


Tiago Peixoto's avatar
Tiago Peixoto committed
378 379
boost::python::object do_move_sweep(GraphInterface& gi, GraphInterface& bgi,
                                    boost::any& emat, boost::any sampler,
380 381 382 383 384 385 386 387 388 389
                                    boost::any cavity_sampler, boost::any omrs,
                                    boost::any omrp, boost::any omrm,
                                    boost::any owr, boost::any ob,
                                    boost::any olabel, vector<int>& vlist,
                                    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,
390 391
                                    size_t nmerges, boost::any omerge_map,
                                    size_t niter,
392 393
                                    partition_stats_t& partition_stats,
                                    bool verbose, rng_t& rng)
394 395 396 397 398 399 400
{
    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;
401
    typedef property_map_type::apply<int32_t,
402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417
                                     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;
418

419 420
    vmap_t merge_map = any_cast<vmap_t>(omerge_map);

Tiago Peixoto's avatar
Tiago Peixoto committed
421 422 423 424
    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,
                        label, vlist, deg_corr, dense, multigraph, beta,
425
                        sequential, parallel, random_move, c, verbose,
426
                        gi.GetMaxEdgeIndex(), nmerges, niter, merge_map,
427
                        partition_stats, rng, S, nmoves, bgi),
Tiago Peixoto's avatar
Tiago Peixoto committed
428
                       mrs, mrp, mrm, wr, b, placeholders::_1,
429
                       std::ref(emat), sampler, cavity_sampler, weighted))();
Tiago Peixoto's avatar
Tiago Peixoto committed
430
    return boost::python::make_tuple(S, nmoves);
431 432
}

433 434 435 436 437
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,
438
                    VertexIndex vertex_index, size_t B, bool weighted, bool empty) const
439 440
    {
        egroups_manage::build(b, oegroups, esrcpos, etgtpos, eweight, g,
441
                              vertex_index, B, weighted, empty);
442 443 444
    }
};

445 446
boost::any do_build_egroups(GraphInterface& gi, GraphInterface& bgi,
                            boost::any ob, boost::any oeweights,
447
                            boost::any oesrcpos, boost::any oetgtpos,
448
                            bool weighted, bool empty)
449 450 451 452 453 454 455
{
    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;
456
    typedef property_map_type::apply<int32_t,
457 458 459 460 461 462 463 464 465
                                     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
466 467 468 469 470
    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()),
471 472
                             placeholders::_1, bgi.GetVertexIndex(),
                             bgi.GetNumberOfVertices(), weighted, empty))();
473 474 475
    return oegroups;
}

476
boost::any do_init_neighbour_sampler(GraphInterface& gi, boost::any oeweights,
477
                                     bool self_loops, bool empty)
478 479 480 481 482 483 484
{
    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
485 486
    run_action<graph_tool::detail::all_graph_views, boost::mpl::true_>()
        (gi, std::bind(init_neighbour_sampler(), placeholders::_1, eweights,
487
                       self_loops, empty, std::ref(osampler)))();
488 489 490
    return osampler;
}

491 492 493 494
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
495
                    std::tuple<boost::any, GraphInterface&> abg) const
496
    {
497
        Graph& bg = *any_cast<Graph*>(get<0>(abg));
498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513
        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
514 515 516 517
    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)))();
518 519
}

Tiago Peixoto's avatar
Tiago Peixoto committed
520 521
boost::python::tuple do_bethe_entropy(GraphInterface& gi, size_t B, boost::any op,
                                      boost::any opv)
522 523 524 525 526 527 528 529 530 531 532
{
    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
533 534 535 536 537
    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);
538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559
}


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
560 561 562
    run_action<graph_tool::detail::all_graph_views, boost::mpl::true_>()
        (gi, std::bind(collect_vertex_marginals_dispatch(),
                       placeholders::_1, b, placeholders::_2),
563 564 565
         vertex_scalar_vector_properties())(op);
}

566 567

struct get_deg_entropy_term
568
{
569 570 571 572 573 574 575 576 577 578 579 580 581 582
    template <class Graph>
    void operator()(Graph& g, double& S) const
    {
        for(auto v : vertices_range(g))
        {
            S -= lgamma_fast(out_degree(v, g) + 1);
            S -= lgamma_fast(in_degreeS()(v, g) + 1);
        }
    }
};

struct get_deg_entropy_term_overlap
{

583
    template <class Graph, class Vprop>
584 585
    void operator()(Graph& g, Vprop b, overlap_stats_t& overlap_stats, size_t N,
                    double& S) const
586 587
    {
#ifdef HAVE_SPARSEHASH
588
        typedef dense_hash_map<int, int, std::hash<int>> map_t;
589
#else
590
        typedef unordered_map<int, int> map_t;
591
#endif
592 593

        map_t in_hist, out_hist;
594 595

#ifdef HAVE_SPARSEHASH
596 597
        in_hist.set_empty_key(numeric_limits<int>::max());
        out_hist.set_empty_key(numeric_limits<int>::max());
598 599
#endif

600
        for (size_t v = 0; v < N; ++v)
601
        {
602 603
            in_hist.clear();
            out_hist.clear();
604

605 606
            auto& half_edges = overlap_stats.get_half_edges(v);
            for (size_t u : half_edges)
607
            {
608 609
                in_hist[b[u]] += in_degreeS()(u, g);
                out_hist[b[u]] += out_degree(u, g);
610
            }
611 612 613 614 615

            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);
616 617 618 619
        }
    }
};

620 621
double do_get_deg_entropy_term(GraphInterface& gi, boost::any ob,
                               overlap_stats_t& overlap_stats, size_t N)
622 623 624 625 626 627
{
    typedef property_map_type::apply<int32_t,
                                     GraphInterface::vertex_index_map_t>::type
        vmap_t;

    double S = 0;
628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643

    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
    {
        run_action<>()
            (gi, std::bind(get_deg_entropy_term(),
                           placeholders::_1, std::ref(S)))();
    }

644 645 646 647
    return S;
}


648 649 650 651 652
vector<int32_t> get_vector(size_t n)
{
    return vector<int32_t>(n);
}

653
template <class Value>
Tiago Peixoto's avatar
Tiago Peixoto committed
654
void vector_map(boost::python::object ovals, boost::python::object omap)
655 656
{

657 658
    multi_array_ref<Value,1> vals = get_array<Value,1>(ovals);
    multi_array_ref<Value,1> map = get_array<Value,1>(omap);
659 660 661 662

    size_t pos = 0;
    for (size_t i = 0; i < vals.size(); ++i)
    {
663
        Value v = vals[i];
664 665 666 667 668 669
        if (map[v] == -1)
            map[v] = pos++;
        vals[i] = map[v];
    }
}

670
template <class Value>
671 672 673
void vector_continuous_map(boost::python::object ovals)
{

674 675
    multi_array_ref<Value,1> vals = get_array<Value,1>(ovals);
    unordered_map<Value, size_t> map;
676 677 678

    for (size_t i = 0; i < vals.size(); ++i)
    {
679
        Value v = vals[i];
680 681 682 683 684 685 686
        auto iter = map.find(v);
        if (iter == map.end())
            iter = map.insert(make_pair(v, map.size())).first;
        vals[i] = iter->second;
    }
}

687
template <class Value>
Tiago Peixoto's avatar
Tiago Peixoto committed
688
void vector_rmap(boost::python::object ovals, boost::python::object omap)
689 690
{

691 692
    multi_array_ref<Value,1> vals = get_array<Value,1>(ovals);
    multi_array_ref<Value,1> map = get_array<Value,1>(omap);
693 694 695 696 697 698 699

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

700 701 702 703 704 705 706 707 708 709 710
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
{
    template <class Graph, class Vprop, class Eprop>
    void operator()(Graph& g, Vprop b, Eprop eweight, size_t N, size_t B,
711
                    bool edges_dl, partition_stats_t& partition_stats) const
712
    {
713
        partition_stats = partition_stats_t(g, b, eweight, N, B, edges_dl);
714 715 716 717 718
    }
};

partition_stats_t
do_get_partition_stats(GraphInterface& gi, boost::any ob, boost::any aeweight,
719
                       size_t N, size_t B, bool edges_dl)
720 721 722 723 724 725 726 727 728 729 730 731 732 733
{
    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;

    partition_stats_t partition_stats;

    vmap_t b = any_cast<vmap_t>(ob);
    emap_t eweight = any_cast<emap_t>(aeweight);

    run_action<>()(gi, std::bind(get_partition_stats(),
734
                                 placeholders::_1, b, eweight, N, B, edges_dl,
735 736 737 738
                                 std::ref(partition_stats)))();
    return partition_stats;
}

739 740 741 742
void export_blockmodel()
{
    using namespace boost::python;

743 744 745 746 747 748 749
    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);

750 751 752 753 754 755
    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);
756 757 758 759 760
    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>);

761
    def("get_vector", get_vector);
762 763 764 765 766 767
    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>);
768 769

    def("create_emat", do_create_emat);
770
    def("create_ehash", do_create_ehash);
771
    def("build_egroups", do_build_egroups);
772
    def("init_neighbour_sampler", do_init_neighbour_sampler);
773 774 775 776

    def("move_sweep", do_move_sweep);

    def("entropy", do_get_ent);
777
    def("entropy_dense", do_get_ent_dense);
778 779
    def("entropy_parallel", do_get_ent_parallel);
    def("deg_entropy_term", do_get_deg_entropy_term);
780 781 782 783 784

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