multiflip: remove O(N^2) corner case with coalescence

src/graph/inference/blockmodel/graph_blockmodel.hh

@@ -128,7 +128,7 @@ public:
           _vweight(uncheck(__avweight, typename std::add_pointer<vweight_t>::type())),
           _eweight(uncheck(__aeweight, typename std::add_pointer<eweight_t>::type())),
           _degs(uncheck(__adegs, typename std::add_pointer<degs_t>::type())),
-          _emat(_bg, rng),
+          _emat(_g, _bg, rng),
           _neighbor_sampler(_g, _eweight),
           _m_entries(num_vertices(_bg))
     {
@@ -1104,25 +1104,29 @@ public:
             d.emplace_back(get<0>(kn.first), get<1>(kn.first), kn.second);
     }
 
-    size_t add_block()
-    {
-        size_t r = boost::add_vertex(_bg);
-        _wr.resize(num_vertices(_bg));
-        _mrm.resize(num_vertices(_bg));
-        _mrp.resize(num_vertices(_bg));
-        _wr[r] = _mrm[r] = _mrp[r] = 0;
-        _bclabel.resize(num_vertices(_bg));
-        _brecsum.resize(num_vertices(_bg));
-        _empty_pos.resize(num_vertices(_bg));
-        _candidate_pos.resize(num_vertices(_bg));
-        add_element(_empty_blocks, _empty_pos, r);
-        for (auto& p : _partition_stats)
-            p.add_block();
-        if (!_egroups.empty())
-            _egroups.add_block();
-        if (_coupled_state != nullptr)
-            _coupled_state->coupled_resize_vertex(r);
-        sync_emat();
+    size_t add_block(size_t n = 1)
+    {
+        _wr.resize(num_vertices(_bg) + n);
+        _mrm.resize(num_vertices(_bg) + n);
+        _mrp.resize(num_vertices(_bg) + n);
+        _bclabel.resize(num_vertices(_bg) + n);
+        _brecsum.resize(num_vertices(_bg) + n);
+        _empty_pos.resize(num_vertices(_bg) + n);
+        _candidate_pos.resize(num_vertices(_bg) + n);
+        size_t r = null_group;
+        for (size_t i = 0; i < n; ++i)
+        {
+            r = boost::add_vertex(_bg);
+            _wr[r] = _mrm[r] = _mrp[r] = 0;
+            add_element(_empty_blocks, _empty_pos, r);
+            for (auto& p : _partition_stats)
+                p.add_block();
+            if (!_egroups.empty())
+                _egroups.add_block();
+            if (_coupled_state != nullptr)
+                _coupled_state->coupled_resize_vertex(r);
+        }
+        _emat.add_block(_bg);
         return r;
     }
 

src/graph/inference/blockmodel/graph_blockmodel_emat.hh

@@ -31,8 +31,8 @@ template <class BGraph>
 class EMat
 {
 public:
-    template <class RNG>
-    EMat(BGraph& bg, RNG&)
+    template <class Graph, class RNG>
+    EMat(Graph&, BGraph& bg, RNG&)
     {
         sync(bg);
     }
@@ -52,6 +52,11 @@ public:
         }
     }
 
+    void add_block(BGraph& bg)
+    {
+        sync(bg);
+    }
+
     typedef typename graph_traits<BGraph>::vertex_descriptor vertex_t;
     typedef typename graph_traits<BGraph>::edge_descriptor edge_t;
 
@@ -125,9 +130,9 @@ class EHash
 {
 public:
 
-    template <class RNG>
-    EHash(BGraph& bg, RNG& rng)
-        : _hash_function(num_vertices(bg), _index, rng),
+    template <class Graph, class RNG>
+    EHash(Graph& g, BGraph& bg, RNG& rng)
+        : _hash_function(num_vertices(g), _index, rng),
           _hash(num_vertices(bg), ehash_t(0, _hash_function))
     {
         sync(bg);
@@ -161,6 +166,15 @@ public:
         }
     }
 
+    void add_block(BGraph& bg)
+    {
+        size_t N = _hash.size();
+        _hash_function.resize(num_vertices(bg));
+        _hash.resize(num_vertices(bg), ehash_t(0, _hash_function));
+        for (size_t i = N; i < _hash.size(); ++i)
+            _hash[i].max_load_factor(.3);
+    }
+
     typedef typename graph_traits<BGraph>::vertex_descriptor vertex_t;
     typedef typename graph_traits<BGraph>::edge_descriptor edge_t;
 

src/graph/inference/blockmodel/graph_blockmodel_multiflip_mcmc.hh

@@ -229,7 +229,6 @@ struct MCMC
             _nmoves++;
         }
 
-
         template <class RNG>
         std::tuple<double, double>
         gibbs_sweep(std::vector<size_t>& vs, size_t r, size_t s,
@@ -422,6 +421,13 @@ struct MCMC
 
             size_t pos = 0;
             std::array<size_t, 2> except = {r, s};
+
+            size_t nB = _groups[r].size();
+            if constexpr (!forward)
+                nB += _groups[s].size();
+            if (_state._empty_blocks.size() < nB)
+                _state.add_block(nB - _state._empty_blocks.size());
+
             for (auto v : _groups[r])
             {
                 size_t t;

src/graph/inference/blockmodel/graph_blockmodel_util.hh

@@ -44,7 +44,7 @@ public:
     virtual size_t sample_block(size_t v, double c, double d, rng_t& rng) = 0;
     virtual double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
                                  bool reverse) = 0;
-    virtual size_t add_block() = 0;
+    virtual size_t add_block(size_t n = 1) = 0;
     virtual void add_edge(const GraphInterface::edge_t& e) = 0;
     virtual void remove_edge(const GraphInterface::edge_t& e) = 0;
     virtual void add_edge_rec(const GraphInterface::edge_t& e) = 0;

src/graph/inference/layers/graph_blockmodel_layers.hh

@@ -957,9 +957,9 @@ struct Layers
             BaseState::set_vertex_weight(v, w);
         }
 
-        size_t add_block()
+        size_t add_block(size_t n = 1)
         {
-            return BaseState::add_block();
+            return BaseState::add_block(n);
             // for (size_t l = 0; l < _layers.size(); ++l)
             // {
             //     auto& state = _layers[l];

src/graph/inference/overlap/graph_blockmodel_overlap.hh

@@ -87,7 +87,7 @@ public:
         : OverlapBlockStateVirtualBase<Ts...>(std::forward<ATs>(args)...),
           _bg(boost::any_cast<std::reference_wrapper<bg_t>>(__abg)),
           _c_mrs(_mrs.get_checked()),
-          _emat(_bg, rng),
+          _emat(_g, _bg, rng),
           _egroups_update(true),
           _overlap_stats(_g, _b, _half_edges, _node_index, num_vertices(_bg)),
           _coupled_state(nullptr)
@@ -1137,25 +1137,29 @@ public:
         }
     }
 
-    size_t add_block()
-    {
-        size_t r = boost::add_vertex(_bg);
-        _wr.resize(num_vertices(_bg));
-        _mrm.resize(num_vertices(_bg));
-        _mrp.resize(num_vertices(_bg));
-        _wr[r] = _mrm[r] = _mrp[r] = 0;
-        _bclabel.resize(num_vertices(_bg));
-        _empty_pos.resize(num_vertices(_bg));
-        _candidate_pos.resize(num_vertices(_bg));
-        add_element(_empty_blocks, _empty_pos, r);
-        _overlap_stats.add_block();
-        for (auto& p : _partition_stats)
-            p.add_block();
-        if (!_egroups.empty())
-            _egroups.add_block();
-        if (_coupled_state != nullptr)
-            _coupled_state->coupled_resize_vertex(r);
-        sync_emat();
+    size_t add_block(size_t n = 1)
+    {
+        _wr.resize(num_vertices(_bg) + n);
+        _mrm.resize(num_vertices(_bg) + n);
+        _mrp.resize(num_vertices(_bg) + n);
+        _bclabel.resize(num_vertices(_bg) + n);
+        _empty_pos.resize(num_vertices(_bg) + n);
+        _candidate_pos.resize(num_vertices(_bg) + n);
+        size_t r = null_group;
+        for (size_t i = 0; i < n; ++i)
+        {
+            r = boost::add_vertex(_bg);
+            _wr[r] = _mrm[r] = _mrp[r] = 0;
+            add_element(_empty_blocks, _empty_pos, r);
+            _overlap_stats.add_block();
+            for (auto& p : _partition_stats)
+                p.add_block();
+            if (!_egroups.empty())
+                _egroups.add_block();
+            if (_coupled_state != nullptr)
+                _coupled_state->coupled_resize_vertex(r);
+        }
+        _emat.add_block(_bg);
         return r;
     }
 

src/graph/inference/partition_centroid/graph_partition_centroid.hh

@@ -161,9 +161,8 @@ public:
         return _wr[_b[v]] - 1;
     }
 
-    constexpr size_t add_block()
+    constexpr void add_block(size_t)
     {
-        return 0;
     }
 
     double virtual_move(size_t v, size_t r, size_t nr)

src/graph/inference/partition_modes/graph_partition_mode_clustering.hh

@@ -259,9 +259,8 @@ public:
         return _wr[_b[v]] - 1;
     }
 
-    constexpr size_t add_block()
+    constexpr void add_block(size_t)
     {
-        return 0;
     }
 
     double virtual_move(size_t v, size_t r, size_t nr)

src/graph_tool/inference/blockmodel.py

@@ -1495,7 +1495,7 @@ class BlockState(object):
                                                           _get_rng())
 
 
-    def multiflip_mcmc_sweep(self, beta=1., c=1., psingle=100, psplit=1,
+    def multiflip_mcmc_sweep(self, beta=1., c=1., psingle=None, psplit=1,
                              pmerge=1, pmergesplit=1, d=0.01, gibbs_sweeps=10,
                              niter=1, entropy_args={}, accept_stats=None,
                              verbose=False, **kwargs):
@@ -1513,8 +1513,9 @@ class BlockState(object):
             node and their block connections; for :math:`c\to\infty` the blocks
             are sampled randomly. Note that only for :math:`c > 0` the MCMC is
             guaranteed to be ergodic.
-        psingle : ``float`` (optional, default: ``100``)
-            Relative probability of proposing a single node move.
+        psingle : ``float`` (optional, default: ``None``)
+            Relative probability of proposing a single node move. If ``None``,
+            it will be selected as the number of nodes in the graph.
         psplit : ``float`` (optional, default: ``1``)
             Relative probability of proposing a group split.
         pmergesplit : ``float`` (optional, default: ``1``)
@@ -1550,6 +1551,8 @@ class BlockState(object):
 
         """
 
+        if psingle is None:
+            psingle = self.g.num_vertices()
         gibbs_sweeps = max(gibbs_sweeps, 1)
         nproposal = Vector_size_t(4)
         nacceptance = Vector_size_t(4)