Detecting mesoscale structures by surprise

The importance of identifying the presence of mesoscale structures in complex networks can be hardly overestimated. So far, much attention has been devoted to the detection of communities, bipartite and core-periphery structures on binary networks: such an effort has led to the definition of a unified framework based upon the score function called surprise, i.e. a p-value that can be assigned to any given partition of nodes, on both undirected and directed networks. Here, we aim at making a step further, by extending the entire framework to the weighted case: after reviewing the application of the surprise-based formalism to the detection of binary mesoscale structures, we present a suitable generalization of it for detecting weighted mesoscale structures, a topic that has received much less attention. To this aim, we analyze four variants of the surprise; from a technical point of view, this amounts at employing four variants of the hypergeometric distribution: the binomial one for the detection of binary communities, the multinomial one for the detection of binary"bimodular"structures and their negative counterparts for the detection of communities and"bimodular"structures on weighted networks. On top of that, we define two"enhanced"variants of surprise, able to encode both binary and weighted constraints and whose definition rests upon two suitable generalizations of the hypergeometric distribution itself. As a result, we present a general, statistically-grounded approach to detect mesoscale structures on networks via a unified, suprise-based framework. To illustrate the performance of our methods, we, first, test them on a variety of well-established, synthetic benchmarks and, then, apply them to several real-world networks, i.e. social, economic, financial and ecological ones. Moreover, we attach to the paper a Python code implementing all the considered variants of surprise.


INTRODUCTION
The importance of identifying the signature of some kind of mesoscopic organization in complex networks (be it due to the presence of communities or bipartite, coreperiphery, bow-tie structures) can be hardly overestimated [1,2], the best example of complex systems whose behavior is deeply affected by their mesoscopic structural organization (e.g. resilience to the propagation of shocks, to the failure of nodes, etc.) being provided by financial networks [3][4][5][6].
So far, much attention has been devoted to the detection of binary mesoscale structures, i.e. communities and, to a far less extent, core-periphery structures: the efforts to solve these problems have led to a number of approaches that are briefly sketched below (for a detailed review of them, see [7,8]).
Community detection has been initially approached by attempting a definition of 'communities' based on the concepts of clustering, cliques, k-core; core-periphery structures have, instead, been defined in a purely topdown fashion, by imagining a fully connected subgraph (i.e. the core) surrounded by (peripherical) vertices ex-clusively linked to the first ones [3]. As stressed in [5], the deterministic character of these definitions makes their tout court application to real-world systems extremely difficult. This is the reason why the intuitive requirements that 'the number of internal edges is larger than the number of external edges' and that 'the core portion of a network is densely connected, while its periphery is loosely connected' [8] are, now, interpreted in a purely probabilistic way: a community has, thus, become a subgraph whose vertices have a larger probability to be interconnected than to be connected to any other vertex in the graph -and analogously for the core-periphery structure. In other words, the top-down approach defining a golden standard and looking for (deviations from) it has left the place to a bottom-up approach where structures are supposed to emerge as the result of non-trivial (i.e. non-casual) interactions between the nodes. cedure for selecting the best model among the ones providing competing descriptions of the data. This has led to the identification of three broad classes of algorithms 1 : according to our intuition, all of them implement some kind of statistical inference, the major difference lying in the way the corresponding test of hypothesis is implemented; from a practical point of view, instead, all these methods are designed for optimization, the functional form of the specific score function determining the class to which a given algorithm belongs.
The most representative algorithms among those belonging to the first class are the ones based on modularity. Although these methods compare the empirical network with a benchmark, they do not provide any indication of the statistical significance of the recovered partition, the reason being that none of them is designed as a proper statistical test. For instance, let us consider the definition of modularity, whose generic addendum is proportional to the term (a ij − p ij ): while it embodies a comparison between the (empirical) adjacency matrix A and the matrix of probability coefficients P defining the benchmark, it does not implement any proper test of hypothesis. The algorithms prescribing to maximize a plain likelihood function belong to this group as well. Although popular, this way of proceeding is known to be affected by overfitting issues: an example is provided by the straight maximization of the likelihood defining the Stochastic Block Model (SBM), or its degree-corrected version, over the whole set of possible partitions, outputting the trivial one where each vertex is a cluster on its own [11].
The aforementioned, major limitation is overcome by the algorithms belonging to the second class. They implement tests of hypothesis either 'a la Fisher orà la Neyman-Pearson, i.e. either defining a single benchmark (the null hypothesis of the first scenario) or two, alternative ones (the null and the alternative hypothesis of the second scenario): from a practical point of view, such a result is achieved by identifying the aforementioned benchmarks with proper probability distributions and the best partition of nodes with the one minimizing the corresponding p-value. Surprise-based algorithms belong to this second group: as it has been shown in [12] for a particular case; such a result will be generalized in what follows -optimizing (asymptotic) surprise amounts at carrying out a (sort of) Likelihood Ratio Test aimed at choosing between two alternative models.
Hypothesis testing can be further refined by allowing for more than two hypotheses to be tested at a time: results of the kind are particularly useful for model selection and, in fact, have produced a plethora of criteria (e.g. the Akaike Information Criterion, the Bayesian Information Criterion and the Minimum Description Length) for singling out the best statistical model out of a basket of competing ones. Generally speaking, optimization, here, seeks for the maximum of a 'corrected' likelihood function embodying the trade-off between accuracy and parsimony of a description. An example of the algorithms belonging to this third class is represented by Infomap; other examples are represented by the recipes employing the SBM within a Bayesian framework (see [13] and the references therein).
With this paper, we pose ourselves within the second research line and adopt a bottom-up approach that prescribes to compare any empirical network structure with the outcome of a properly-defined benchmark model. To this aim, we devise a unified framework for mesoscale structures detection based upon the score function called surprise 2 , i.e. a p-value that can be assigned to any given partition of nodes, on both undirected and directed networks: while for binary community detection this is achieved by employing the binomial hypergeometric distribution [15,39], with 'bimodular' structures like the bipartite and the core-periphery ones, one needs to consider its multinomial variant [12]. Here, however, we aim at making a step further, by extending the entire framework to the weighted case. As a result, we present a general, statistically-grounded approach to the problem of detecting mesoscale structures on networks via a unified, suprise-based framework.
Before moving forward, we would like to stress that the use of the hypergeometric distribution to carry out tests of hypothesis on networks is not novel: examples are provided by the papers [16] (where the authors introduce a method to provide a statistically-validated, monopartite projection of a bipartite network -the considered null hypothesis encoding the heterogeneity of the system), [17] (where the authors employ the same validation procedure to detect cores of communities within each set of nodes of a bipartite system) and [18] (where the authors extend the framework proposed in the aforementioned references to carry out a statistical validation of motifs observed in hypergraphs). For a review on the use of the hypergeometric distribution for network analyses see [19] and the references therein.

MESOSCALE STRUCTURES DETECTION VIA EXACT TESTS
Surprise has recently received a lot of attention: the advantages of employing such a score function have been extensively discussed in [12,15,[20][21][22][23]39] where researchers have tested and compared its performance from a purely numerical perspective. However, a characterization of the statistical properties of surprise is still missing: for this reason, we will, first, make an effort to 'translate' the problem of detecting a given mesoscale network structure into a proper exact significance test (to the best of our knowledge, the only other attempt of the kindspecifically, to detect communities -is the one in [24]) and, then, show how the rich -yet, still underexplored -surprise-based formalism can properly answer such a question.
The basic equation underlying exact tests reads (1) and returns the probability of observing an outcome of the random variable X which is 'more extreme' than the realized one, i.e. x * . In the setting above, f represents the distribution encoding the null hypothesis and Pr(x ≥ x * ) -commonly known with the name of p-value -answers the question is the realized value X = x * compatible with the (null) hypothesis that X is distributed according to f ? The p-value constitutes the basic quantity for carrying out any significance test: hence, in what follows we will tackle the problem of detecting the signature of statistically-significant mesoscale organizations by individuating specific tests (or, equivalently, suitable functional forms for f ).

DETECTION OF MESOSCALE STRUCTURES IN BINARY NETWORKS
Community detection. The topic of nodes partitioning into densely connected groups has traditionally received a lot of attention, with applications ranging from the analysis of social networks to the definition of novel recommendation systems [1,2]. Within the surprisebased framework, binary community detection is carried out via the identification i.e. by calculating the p-value of a binomial hypergeometric distribution whose parameters read as above. In this formalism, the • subscript will be meant to indicate quantities that are internal to communities, while the • subscript will be meant to indicate quantities that are external to communities. More precisely, the binomial coefficient V• l• enumerates the number of ways l • links can be redistributed within communities, i.e. over the available V • node pairs, while the binomial coefficient V• l• enumerates the number of ways the remaining l • = L − l • links can be redistributed between communities, i.e. over the remaining V • = V − V • node pairs. Notice that, although l • is 'naturally' bounded by the value V • , it cannot exceed L -whence the usual . From a merely statistical point of view, surprise 'considers' a network as a population of V node pairs, L of which have been drawn; out of the L extracted ones, l • node pairs have the desired feature of being 'internal' to communities, since they connect some of the node pairs belonging to the set V • . Hence, for a given partition of nodes into communities, S quantifies the probability of observing at least l * • 'successes' (i.e. intra-cluster edges) out of L draws: the lower this probability, the 'more surprising' the observation of the corresponding partition, hence the 'better' the partition itself.
Asymptotic results. We can gain more insight into the surprise-based formalism above upon deriving an asymptotic expression for S [25]. To this aim, let us consider that it can be simplified upon Stirling-approximating the binomial coefficients that appear within it. By exploiting the recipe n! √ 2πn n e n , S can be rewritten as where the expression defines a Bernoulli probability mass function, the parameters appearing in eq. (4) read p = L V and p i = li Vi and the coefficient in front of the sum is . Equation (4) makes it explicit that employing S for binary community detection ultimately amounts at comparing the description of a networked configuration provided by the Random Graph Model (RGM), and encoded into the expression with the description of the same configuration provided by the Stochastic Block Model (SBM) [26], and encoded into the expression Naturally, the SBM takes as input the two sets indexed by • and • and distinguish the connections found 'within the clusters' -contributing to the probability of the whole configuration with the term from the ones found 'between the clusters' -contributing to the probability of the whole configuration with the term Ber Notice also that the asymptotic expression of the surprise guarantees that the parameters of the null models defining it are tuned according to the maximum-of-thelikelihood principle. To see this explicitly, let us consider Ber(x, y, z) whose log-likelihood reads upon maximizing it with respect to z, one finds z = y x . The way S works, i.e. by comparing two different null models, is reminiscent of more traditional likelihood ratio tests, where a null hypothesis H 0 (in our case, a given partition is compatible with the RGM) is tested against an alternative hypothesis H 1 (in our case, the given partition is compatible with the SBM): as the asymptotic expression of the surprise clarifies, minimizing it amounts at finding the partition least likely to occour under the RGM than under the SBM.
'Bimodular' structures detection. The surprisebased framework can be easily extended to detect what can be called 'bimodular structures', a term that will be used to compactly indicate core-periphery [27][28][29] and bipartite structures [30,31]. The reason for adopting such a terminology lies in the evidence that both kinds of structures are defined by bimodular partitions, i.e. partitions of nodes into two different groups.
As shown elsewhere [12], the issue of detecting binary 'bimodular' structures can be addressed by considering a multivariate (or multinomial ) hypergeometric distribution, i.e. by identifying indicates the number of node pairs between the modules • and • and l ≡ L − (l • + l • ) indicates the number of links that must be assigned therein. While the binomial coefficient V• l• enumerates the number of ways l • links can redistributed within the first module (e.g. the core portion) and the binomial coefficient V• l• enumerates the number of ways l • links can redistributed within the second module (e.g. the periphery portion), the third binomial enumerates the number of ways the remaining L−(l • +l • ) links can be redistributed between the first and the second module, i.e. over the remaining This choice induces the definition of the binary bimodular surprise analogously to the univariate case, l • and l • are 'naturally' bounded by the values V • and V • -notice, however, that the sum l • + l • cannot exceed L (although it may not reach such a value, e.g. in case V • + V • < L).
Asymptotic results. Analogously to the univariate case, the asymptotic expression for S can be derived upon Stirling-approximating the binomial coefficients appearing within it. This time, the recipe n! √ 2πn n e n leads to where, as before, Ber(x, y, z) = z y (1 − z) x−y defines a Bernoulli probability mass function and the parameters read p = L V and p i = li Vi ; the numerical coefficient appearing in front of the whole expression, now, reads The quantity S compares the description of a networked configuration provided by the RGM, and encoded into the expression Ber(V, L, p) = p L (1 − p) V −L , with the description of the same configuration provided by the SBM (now, defined by three -instead of two -different blocks), 'represented' by the denominator of the expression defined in eq. (12), i.e. i=•,•,

DETECTION OF MESOSCALE STRUCTURES IN WEIGHTED NETWORKS
Community detection. Within the surprise-based framework, the problem of detecting binary communities has been rephrased as an aleatory experiment whose random variable is the number of links within communities. Interestingly enough, such an experiment can be easily mapped into a counting problem, allowing us to interpret S as indicating the number of configurations whose number of 'internal' links (i.e. within communities) is larger than the observed one. When dealing with weighted networks, we would like to proceed along similar guidelines and consider the total, 'internal' weight as our new random variable, to be redistributed across the available node pairs. Adopting this approach has four major consequences: 1) weights must be considered as composed by an integer number of binary links, 2) each node pair must be allowed to be occupied by more than one link and 3) the total weight must be allowed to vary even beyond the network size (when handling real-world networks, the case W V is often encountered).
In this case, the proper setting to define an aleatory experiment satisfying the requests above is provided by the so-called stars and bars model, a combinatorial technique that has been introduced to handle the counting of configurations with multiple occupancies. Basically, the problem of counting in how many ways w • particles (our links) can be redistributed among V • boxes (our node pairs), while allowing more than one particle to occupy each box, can be tackled by allowing both the particles and the bars 'delimiting' the boxes to be permuted [32]. Since V • boxes are delimited by V • − 1 bars, a term like is needed. In order to better grasp the meaning of such a term, let us make a simple example. Let us imagine to observe a network with three nodes and two links, carrying a weight of 1 and 2, respectively. Now, were we interested in a purely binary analysis, we may ask ourselves in how many ways we could place the two links among the N (N −1) 2 = 3(3−1) 2 = 3 available pairs: the answer is provided by the 'binary' binomial coefficient V• l• = 3 2 = 3. The implicit assumption we make is that the three links must not occupy the same node pairs -otherwise the total number of connections wouldn't be preserved.
This perspective changes from the purely weighted point of view. Since we are now interested in preserving just the total weight of our network, irrespectively of the number of connections it is placed upon, the number of admissible configurations amounts precisely at . Such a number is larger than before since, now, weights are 'disaggregated' into binary links and multiple occupations of the latter ones are allowed.
The considerations above lead us to generalize the community detection problem to the weighted case by identifying i.e. by replacing the binomial hypergeometric distribution considered in the purely binary case with a negative hypergeometric distribution, a choice inducing the definition of the weighted surprise where the binomial coefficient V enumerates the number of ways w • links can be redistributed within communities, i.e. over the available V • node pairs, and the binomial coefficient enumerates the number of ways the remaining w • = W − w • links can be redistributed between communities, i.e. over the remaining V • = V − V • node pairs. Differently from the binary case, the sum ranges up to the maximum empirical weight of the network, i.e. w • ∈ [w * • , W ].
Asymptotic results. The asymptotic expression for W can be deduced by following the same reasoning that has allowed us to derive the asymptotic expression for S . Stirling-approximating the binomial coefficients entering into the definition of W leads to the writing where the expression defines a geometric probability mass function and the parameters appearing in eq. (17) read q = W V +W −1 and q i = wi Vi+wi−1 . In the weighted case, the Bernoulli probability mass function appearing in the asymptotic expression of S is replaced by a geometric probability mass function: this implies that (asymptotically) the comparison is, now, carried out between the description of a networked configuration provided by the Weighted Random Graph Model (WRGM), and encoded into the expression with the description of the same configuration provided by the Weighted Stochastic Block Model (WSBM) and encoded into the expression As for the binary case, the asymptotic expression of the weighted surprise clarifies that minimizing it amounts at finding the partition least likely to occour under the WRGM than under the WSBM.
As in the binary case, the parameters characterizing the geometric probability mass functions defining the asymptotic weighted surprise can be estimated via the maximum-of-the-likelihood principle, according to which the log-likelihood of the expression Geo(x, y, z) reads upon maximizing it with respect to z, one finds z = y x+y : hence, one can pose q W V +W and q i wi Vi+wi . The numerical coefficient appearing in front of the whole expression, instead, reads with µ V q and µ i V i q i .
'Bimodular' structures detection. Let us now introduce the third generalization of the suprise-based formalism: following the same line of reasoning that led us to approach the detection of binary 'bimodular' structures by considering the multinomial analogoue of the distribution introduced for binary community detection, we focus on the multinomial (or multivariate) negative hypergeometric distribution, i.e.
while the binomial coefficient V•+w•−1 w• enumerates the number of ways w • links can redistributed within the first module (e.g. the core portion), the binomial coefficient V•+w•−1 w• enumerates the number of ways w • links can be redistributed within the second module (e.g. the periphery portion) and the binomial coeffi- enumerates the number of ways the remaining w ≡ W − (w • + w • ) links can be redistributed between the first and the second module, i.e. over the remaining V ≡ V − (V • + V • ) node pairs. Such a position induces the definition of the weighted bimodular surprise as for the (weighted) community detection, weights are understood as integer numbers -equivalently, as composed by an integer number of binary links. For what concerns the limits of the summations, w • and w • are 'naturally' bounded by W ; notice, however, that the sum w • + w • itself cannot exceed such a value.
Asymptotic results. Let us now derive the asymptotic expression for W . As usual, let us Stirling-approximate the binomial coefficients entering into the definition of W ; such a simplification leads us to the expression where, as before, Geo(x, y, z) = z y (1 − z) x defines a geometric probability mass function and the parameters read q = W V +W −1 and q i = wi Vi+wi−1 but can be approximated, according to the maximum-of-the-likelihood principle, as q W V +W and q i wi Vi+wi ; the numerical coefficient multiplying the whole expression reads with µ V q and µ i V i q i . The analysis of W in the asymptotic regime reveals that it compares the description of a networked configuration provided by the WRGM, and encoded into the expression Geo(V, W, q) = q W (1 − q) V , with the description of the same configuration provided by the WSBM (now, defined by three -instead of two -different blocks), 'represented' by the denominator of the expression defined in eq. (25), i.e. i=•,•,

ENHANCED DETECTION OF MESOSCALE STRUCTURES IN NETWORKS
Community detection. The recipe to detect communities on weighted networks can be further refined to account for the information encoded into the total number of links, beside the one provided by the total weight. Generally speaking, this can be realized by 'combining' two of the distributions introduced above. To this aim, let us proceed in a two-step fashion: first, let us recall that the number of ways L links can be placed among V node pairs, in such a way that l • connections are 'internal' to the clusters while the remaining L − l • ones are, instead, 'external' is precisely now, for each of the binary configurations listed above, W − L links remain to be assigned: while w • − l • of them must be placed within the clusters, on top of the l • available 'internal' links, the remaining (W − L) − (w • − l • ) ones must be placed between the clusters, on top of the L−l • available inter-cluster connections. Hence, the 'conditional' negative hypergeometric distribution reading remains naturally defined; now, 'combining' the two distributions above, simplifying and re-arranging, the generic term of the enhanced hypergeometric distribution can be rewritten as with a clear meaning of the symbols. An analytical characterization of it is provided into the Appendices: for the moment, let us simply notice that the definition provided above works for the values 0 < l • < L. By posing f (l • , w • ) ≡ EH(l • , w • |V, V • , L, W ), our novel distribution induces the definition of the enhanced surprise, i.e.
although l • and w • − l • are 'naturally' bounded by V • and W −L, respectively, the former one cannot exceed L.
In order to better understand how the enhanced surprise works, let us consider again the aforementioned example: given a network with three nodes and two links, carrying a weight of 1 and 2, respectively, we observe V• l• Let us imagine to observe a network with three nodes and two links, carrying a weight of 1 and 2, respectively. Were we interested in a purely binary analysis, we may ask ourselves in how many ways we could place the two links among the 3 available pairs: the answer is provided by the 'binary' binomial coefficient V• l• = 3 2 = 3. Were we interested in a purely weighted analysis, we may ask ourselves in how many ways we could place the two links among the 3 available pairs while preserving the total weight of our network, irrespectively of the number of connections it is placed upon; the number of admissible configurations becomes V•+w•−1 w• = 2+3 3 = 10. Such a number is larger than before since, now, weights are 'disaggregated' into binary links and multiple occupations of the latter ones are allowed. Were we interested in the enhanced analysis, we may ask ourselves in how many ways we could place the two links among the 3 available pairs while preserving both the total number of links and the total weight of the network; the number of admissible configurations becomes tions with exactly the same total weight. If we, now, constrain both the total number of links and the total weight of the network, the number of admissible config- as it can be easily verified upon explicitly listing them. Naturally, the configurations 'admissible' by the enhanced surprise are a subset of the configurations 'admissible' by the weighted surprise, i.e. precisely the ones with the desired number of links (see also fig. 1).
Asymptotic results. Analogously to the other functionals, the asymptotic expression of E can be derived by Stirling-approximating the binomial coefficients entering into its definition as well: (32) in the formula above, the expression defines a Bose-Fermi probability mass function [33]. While a Bernoulli probability mass function characterizes the asymptotic behavior of the purely binary surprise and a geometric probability mass function characterizes the asymptotic behavior of the purely weighted surprise, the enhanced surprise is asymptotically characterized by a distribution whose functional form is halfway between the two previous ones: while its Bernoulli-like portion controls for the 'presence' of the links, its 'conditional' geometric-like portion controls for the magnitude of the 'remaining' weights. To stress its 'mixed' character, such a distribution has been named 'Bose-Fermi': remarkably, it can be retrieved within the Exponential Random Graphs framework as a consequence of Shannon entropy maximization, constrained to simultaneously reproduce both the total number of links and the total weight of a network [33,34]. The parameters appearing in eq. (32) read p = L V , while the first class of parameters can be tuned according to the maximum-ofthe-likelihood principle, the ones belonging to the second class can be approximated according to the same recipe. To see this explicitly, let us consider BF(x, y, u, z, t), whose log-likelihood reads L = y ln z+(x−y) ln(1−z)+(u−y) ln t+y ln(1−t); (34) upon maximizing it with respect to z, one finds z = y x ; upon maximizing it with respect to t one finds t = u−y u -hence, one can pose r W −L W , r i wi−li wi . The numerical coefficient multiplying the whole expression is defined as Similarly to the other cases, the asymptotic expression of the enhanced suprise compares the description of a networked configuration provided by the Enhanced Random Graph Model (ERGM), and encoded into the expression with the description of the same configuration provided by its block-wise counterpart, i.e. the Enhanced Stochastic Block Model (ESBM) and encoded into the expression 'Bimodular' structures detection. The last generalization of surprise concerns its use for the detection of 'bimodular' structures within the enhanced framework. This amounts at considering the following 'multinomial variant' of the enhanced hypergeometric distribution indicates the number of node pairs between the modules • and • and l ≡ L − (l • + l • ) indicates the number of links that must be assigned therein. An analytical characterization of it is provided into the Appendices: for the moment, let us simply notice that the definition provided above works for the Notice that l • and l • are 'naturally' bounded by V • and V • : still, their sum cannot exceed L; analogously, w • − l • and w • − l • are 'naturally' bounded by W − L: still, their sum itself cannot exceed W − L.
As for its binomial counterpart, the expression of the MEH can be rearranged in a term-by-term fashion, in such a way that the module-specific binomial coefficients can be grouped together. Upon doing so, it becomes clearer that the MEH counts the number of ways w • − l • links can be placed on top of the l • binary links characterizing the connectance of the • module, times the number of ways w • − l • links can be placed on top of the l • binary links characterizing the connectance of the • module, times the number of ways the remaining ) links can be placed on top of the L − (l • + l • ) binary links characterizing the connectance of the third module.
Asymptotic results. The enhanced bimodular surprise admits an asymptotic expression as well, i.e.

Full expression
Community detection 'Bimodular' structures detection  Table illustrating all the generalizations of the surprise-based formalism proposed in the present paper and showing both the full and the asymptotic expression of each probability mass function -with the only exception of the numerical coefficient characterizing each asymptotic expression. To sum up, detecting a weighted mesoscale structure implies considering the negative version of the probability mass function working in the corresponding binary case (e.g. moving from the hypergeometric to the negative hypergeometric one); detecting a 'bimodular' structure, instead, implies considering the multinomial version of the probability mass function working in the corresponding binary case (e.g. moving from the binomial hypergeometric to the multinomial one). Upon employing the multiset notation, a nice formal symmetry can be recovered between the purely binary and the purely weighted cases.
that can be recovered by Stirling-approximating the binomial coefficients entering into the definition of E . While the expression BF denotes a Bose-Fermi probability mass function [33], the parameters appearing in eq. (38) read p = L V , p i = li Vi and r = W −L W −1 , r i = wi−li wi−1 ; as for E , the maximum-of-the-likelihood principle determines (approximates) the first (second) class of parameters.
The numerical coefficient multiplying the whole expression is defined as where Lr and µ i l i r i . As observed for the other functionals, the asymptotic expression of the enhanced bimodular suprise compares the description of a networked configuration provided by the ERGM, and encoded into the expression with the description of the same configuration provided by its block-wise counterpart, i.e. the ESBM (now, defined by three -instead of two -different blocks), 'represented' by the denominator of the expression defined in eq. (38), Table II gathers all the variants of the surprise-based formalism, illustrating both the full and the asymptotic expression for each of them. To sum up, detecting a weighted mesoscale structure implies considering the negative version of the probability mass function working in the corresponding binary case (e.g. moving from the hypergeometric to the negative hypergeometric one); detecting a 'bimodular' structure, instead, implies considering the multinomial version of the probability mass function working in the corresponding binary case (e.g. moving from the binomial hypergeometric to the multinomial one).

RESULTS
The previous sections have been devoted to the description of the surprise-based formalism for detecting a number of mesoscale structures; let us now test it on a bunch of synthetic and real-world configurations.
Consistency checks. Let us start by checking the consistency of our surprise-based formalism. Since we have rephrased the problem of detecting any mesoscale structure into an exact significance test, the limitations of our formalism are the same ones affecting the tests of the kind. More precisely, any significance test is characterized by a parameter known as type I error rate and usually denoted with α: it quantifies the percentage of times the considered test provides a false positive. In other words, any statistical test is known to 'fail'; in our case, this amounts at recognizing that a significant mesoscale Any exact statistical test is characterized by a parameter, usually denoted with α and known as the type I error rate, that quantifies the percentage of times the test provides a false positive: it can be kept below a given threshold by adjusting the threshold of the p-value of the corresponding test, e.g. by deeming its response as significant in case it is less than 0.05. Over 10.000 networks sampled by the RGM with p = 0.2, a given planted partition -from left to right: two, three and four clusters with different dimensions -is indeed recovered as significant 5% of the times. Notably, since the three upper panels show the CDF of the empirical values of S , they also provide an information about its distribution, which is uniform over the unit interval.
structure can be detected, in case there is none, a percentage α of the times: this number can be kept 'small' by adjusting the threshold of the p-value of the corresponding test -e.g. by deeming its response as significant in case it is less than 0.05.
Remarkably, the aforementioned behavior can be explicitly tested for each of the cases considered above. Whenever exercises like these are carried out, it is of utmost importance to be as clear as possible about the null hypothesis tested: here, we aim at testing whether a given partition is significant or not when the null hypothesis is true. For the sake of illustration, let us consider the problem of detecting communities on binary networks: as we learnt from the asymptotic expression of S , it compares the description of a network provided by the RGM with that of the same network provided by the SBM, thus suggesting the RGM as the model playing the role of H 0 . Hence, let us generate many networks from the RGM and, for each of them, let us impose a partitionthe same for each sampled configuration, over which our SBM is tuned; finally, let us calculate S on each of them.
The results of such an experiment are shown in fig. 2: over 10.000 networks sampled from the RGM, whose only parameter has been set to p = 0.2, the (given) planted partition -i.e. two, three and four clusters with different dimensions -is, indeed, recovered as significant 5% of the times; notice that such a result holds true irrespectively of the details of the partition imposed on the sampled configurations.
We have also repeated such an experiment for the problem of detecting communities on weighted networks: as fig. 3 shows, the same results are recovered, i.e. over 10.000 networks sampled from the WRGM, whose only parameter has been set to q = 0.2, a (given) planted partition of two, three and four clusters with different dimensions is recognized as significant 5% of the times; again, such a result holds true irrespectively of the details of the partition imposed on our sampled configurations. As in the binary case, the type I error rate can be kept below a given threshold by adjusting the threshold of the p-value of the corresponding test, e.g. by deeming its response as significant in case it is less than 0.05. Over 10.000 networks sampled by the WRGM with q = 0.2, a given planted partition -from left to right: two, three and four clusters with different dimensions -is indeed recovered as significant 5% of the times. Notably, since the three upper panels show the CDF of the empirical values of W , they also provide an information about its distribution, which is uniform over the unit interval.
learn that the latter ones are distributed uniformly over the unit interval, i.e. S ∼ U[0, 1] and W ∼ U[0, 1] -an evidence further confirming the behavior of the p-value of any exact significance test.
Checking the behavior of the multivariate versions of our hypergeometric distributions is, instead, much more difficult, the reason lying in the evidence that the null hypothesis is, now, a multivariate one and few results about the behavior of multivariate p-values are known. Still, we can say something about the behavior of the marginal pvalues; for the sake of illustration, let us consider the ones induced by S and W and reading and (43) respectively. As evident from the formulas above, marginal distributions of multivariate hypergeometric distributions are hypergeometric themselves: hence, the behavior of our marginal p-values is expected to match the one observed in the univariate cases, leading to recover and W (•) ∼ U[0, 1]. Analogously, for the enhanced surprise.
Our findings can be also described from a different perspective, i.e. by answering the question what is the expected value of the surprise on a network sampled from the RGM ensemble? To answer this question let us focus on the simplest case of detecting communities on binary networks and consider that, irrespectively from the partition that is planted on the sampled configuration, the For each configuration, 20 cliques have been considered; the size of the latter ones is left to vary as follows: K3, K4, K5, K8, K10, K15, K20. While surprise always recovers the planted partition, modularity misses the partitions with K3 and K4 and Infomap misses the partition with K3 -a result that may be a consequence of the resolution limit which is known to affect the last two algorithms. Bottom row: comparison between the 'ring of binary cliques' and the 'ring of weighted cliques' cases. The result according to which surprise minimization is able to discriminate the cliques linked in a ring-like fashion changes once the weights come into play: in fact, as the weight of the links connecting any two cliques is risen, the algorithm reveals as 'communities' pairs of tightly-connected cliques. expected number of links will be L = pV while the expected number of 'internal' links will be l * • = pV • , where p is the parameter defining the RGM; hence, the surprise becomes (to be fully consistent we should have rounded both the value of L and that of l * • to the nearest integer but our conclusions are still valid). In order to evaluate the expression above, let us consider that the expected value of the hypergeometric distribution defining the surprise reads the validity of the first passage resting upon the evidence that the network is generated via the RGM. Hence, the surprise becomes a sum over the values l • ≥ l • , i.e. those that are more extreme that the average. Since the hypergeometric distribution is peaked around its average, the result above suggests that summing over larger values encodes half of the probability mass, i.e. compatible with those of a uniform distribution whose support coincides with the unit interval.
Comparison with the binary modularity. For the sake of comparison, let us consider the modularity function for the problem of community detection on binary, undirected networks. It is defined as where, for consistency, we employ the RGM as a benchmark, i.e. a ij = p = L V = 2L N (N −1) , ∀ i < j. From the definition of modularity, it follows that where N c and l * c indicate the size of community c and the number of links within it, respectively -naturally, l * • = c l * c and V • = c N c (N c − 1). This calculation clarifies that a positive value of Q implies that the 'internal' probability of connection p • is larger than the one 'predicted' by the RGM, i.e.
while a null modularity value implies that p • matches the only parameter defining the RGM. The result above can be also restated as follows: a null modularity value points out that the mesoscale structure of the configuration at hand has not been generated by a SBM (or, equivalently, does not need a SBM to be explained). Notice, however, that the modularity provides no indication about the statistical significance of the recovered partition.
Comparison with the weighted modularity. A similar conclusion can be reached upon considering the modularity function for community detection on weighted, undirected networks, defined as where we have employed the Weighted Random Graph Model (WRGM) as a benchmark, according to which [35]. From the definition of weighted modularity provided above, it follows that : what differentiates the two pairs of curves is the value of the binary mixing parameter. When focusing on binary, undirected networks, we have considered τ1 = −2 and τ2 = −1 (the average degree is 20 and the maximum degree is 50). When binary, directed configurations are considered, µt refers to in-degrees, which are distributed according to a power-law while the out-degrees are kept constant for all nodes; the other input parameters, instead, are the same used for undirected configurations. When weighted, undirected configurations are considered, an additional mixing parameters is needed, i.e. µw, accounting for the percentage of a node strength to be distributed on the links that connect it to the nodes outside its own community; the exponent of the strength distribution has been set to 1.5 for all realizations considered here. When weighted, directed configurations are considered, µw refers to in-strengths. The trend of the NMI, plotted as a function of the mixing parameters, reveals Infomap to be a strong performer even if its performance decreases abruptly as the mixing parameters exceed a threshold value that depends on the particular setting of the benchmark; the performance of modularity, instead, 'degrades' less rapidly. The performance of surprise constitutes a good compromise between the robustness of modularity and the steadily high accuracy of Infomap.
where N c and w * c indicate the number of nodes of community c and the weight of links within community c, respectively -naturally, w * • = c w * c . As in the binary case, the calculation above clarifies that a positive value of Q implies that the expected weight of any 'internal' link is larger than the one predicted by the WRGM, i.e.
a null modularity value, on the other hand, implies that w ij WSBM coincides with the expectation coming from the WRGM. Analogously to the binary case, the result above can be also restated as follows: a null modularity value points out that the mesoscale structure of the configuration at hand has not been generated by a WSBM (or, equivalently, does not need a WSBM to be explained). Again, however, the modularity provides no indication about the statistical significance of the recovered partition.
Comparing mesoscale structures detection methods. The previous subsections have been devoted to check the consistency of our surprise-based formalism and reveal the 'statistical flaws' of other approaches. Let us now consider two popular algorithms for mesoscale structures detection, i.e. modularity maximization (now, Q has been considered in its 'full' definition, e.g. a ij = p ij = kikj 2L , ∀ i < j for binary, undirected configurations) and Infomap [36], and carry out more systematical comparisons between the former and the surprise. Upon doing so, we are able to compare one algorithm per class, i.e. modularity for the first class, surprise for the second class and Infomap for the third class.
To this aim, we have focused on different kinds of benchmarks, i.e. classes of synthetic networks with welldefined planted partitions, the aim being that of inspecting the goodness of a given algorithm in recovering the imposed partition. As an indicator of the goodness of the partition retrieved by each algorithm, we have followed [17] and employed three different indices. The first one is the normalized mutual information (NMI), defined as where partitions X and Y are compared, H(X) = − x f x ln f x and f x = nx n is the fraction of nodes assigned to the cluster labeled with x; analogously, for the partition Y . The term I(X, Y ) = x y f xy ln fxy fxfy is the 'proper' mutual information and f xy = nxy n is the fraction of nodes assigned to cluster x in partition X and to cluster y in partition Y . Naturally, I(X, Y ) equals 1 if the partitions are identical and 0 if the partitions are independent.
The second index we have considered is the adjusted Rand index (ARI), defined as ARI = T P + T N − T P + T N T P + F P + T N + F N − T P + T N and representing a sort of accuracy 3 'corrected' by a term that quantifies the agreement between the reference partition and a random partition -the term 'random' referring to the Permutation Model; equivalently, the closer the ARI to 0, the more 'random' the provided partition. The third index we have considered is the adjusted Wallace index (AWI), defined as AWI = T P − T P T P + F P − T P and representing a sort of 'corrected' positive predicted value. Again, the closer the AWI to 0, the more 'random' the provided partition.
First, let us inspect the performance of modularity, surprise and Infomap to detect cliques arranged in a ring. Specifically, we have considered 7, different ring-like configurations, each one linking 20 binary cliques (i.e. K 3 , K 4 , K 5 , K 8 , K 10 , K 15 , K 20 ). As fig. 4 reveals, surprise always recovers the planted partition; on the other hand, modularity maximization leads to miss the partitions with K 3 and K 4 and Infomap misses the partition with K 3 , a result that may be a consequence of the resolution limit, affecting both the aforementioned algorithms.
Let us now ask us if the presence of weights affects the detection of mesoscale structures. Generally speaking, the answer is yes, as the comparison between the 'ring of binary cliques' and the 'ring of weighted cliques' cases shows. In particular, the result according to which surprise minimization is able to discriminate the inter-linked cliques changes once the weight of the links connecting any two cliques is risen: in fact, this leads the algorithm to reveal as 'communities' two tightly-connected pairs of cliques, now. We explicitly notice that the results shown in fig. 4 also depend on the relative magnitude of the weights of the intra-cliques and of the inter-cliques links: however, as long as the inter-cliques weight is up to two orders of magnitude larger than the intra-cliques one, it holds true.
In order to expand the set of comparisons, we have focused on two different kinds of well-established benchmarks, i.e. the Lancichinetti-Fortunato-Radicchi (LFR) one and the Aldecoa's 'relaxed-caveman' (RC) one [38? ]. networks, consisting of 512 nodes, grouped in 16 communities arranged in a ring-like fashion and whose sizes obey a power-law whose exponent has been set to 1.8; the smallest community is composed by 3 nodes. Such initial configurations are progressively 'degraded' according to the following mechanism: first, a percentage p of links is selected randomly and removed; afterwards, a percentage p of links is selected randomly and rewired. The trend of the NMI, plotted as a function of the (single) 'degradation' parameter p, driving the evolution of the initial ring-of-cliques towards a progressively less-defined clustered configuration, reveals surprise to outperform both modularity and Infomap. From a more general perspective, these results confirm what has been already observed elsewhere, i.e. that the best-performing algorithms on the LFR benchmarks often perform poorly on the RC benchmarks and vice versa.
The LFR benchmark is a special case of the plantedpartition model, in which groups have different sizes and nodes have different degrees -hence constituting a refinement of the GN benchmark, where groups have equal size and nodes have the same expected degree [38]. When binary, undirected configurations are considered, the degrees of nodes are distributed according to a power-law with exponent τ 1 while the sizes of communities are distributed according to a power-law with exponent τ 2 . Once the sizes of communities have been drawn each node 'receives' its own degree, say k i : the percentage of these links connecting node i with other internal nodes is, then, chosen to be (1−µ t )k i , with µ t playing the role of mixing parameter that controls for the sharpness of the planted partition. For binary, directed configurations, µ t refers to in-degrees, which are distributed according to a powerlaw while the out-degrees are kept constant for all nodes; the other input parameters, instead, are the same used for undirected configurations.
Results on specific implementations of the LFR benchmark are shown in fig. 6. Infomap is, generally speaking, a strong performer; as evident upon looking at the first row, however, its performance decreases abruptly as the mixing parameter exceeds a threshold value that depends on the particular setting of the LFR benchmark. Modularity, instead, seems to be more robust (i.e. its performance 'degrades' less rapidly as µ t increases) although the resolution limit manifests itself when configurations with small communities are considered. Overall, the performance of surprise seems to constitute a good compromise between the robustness of modularity and the steadily high accuracy of Infomap. We would also like to stress that surprise competes with modularity although it employes much less information than the latter: in fact, while the benchmark employed by modularity coincides with the (sparse version of the) Con- . The one we have defined here progressively 'degrades' an initial configuration, defined by 1) a completely connected core, 2) an empty periphery, 3) an intermediate part whose link density amounts at pcp = 0.5. Such a configuration is 'degraded' by progressively filling the periphery and emptying the core. This is achieved by 1) considering all peripherical node pairs and link them with probability q; 2) considering all core node pairs and keep them linked with probability 1 − q: varying q in the interval [0, pcp] allows us to span a range of configurations starting with Borgatti-Everett and ending with Erdös-Rényi. As expected, the performance of the surprise worsens as the 'degradation' parameter becomes closer to pcp = 0.5; however, both the NMI and the ARI indices steadily remain very close to 1. Bottom row: the presence of weights affects the detection of 'bimodular' mesoscale structures as well. In fact, rising the weight of any two links connecting the core with the periphery of a toy network allows the two nodes originally part of the periphery to be detected as belonging to the core; analogously, if a bipartite topology is modified by adding weights between some of the nodes belonging to the same layer, a core-periphery structure will now be detected as significant.
figuration Model -hence, encodes the information on the entire degree sequence -surprise compares the RGM with . Overall, link weights refine the picture provided by just considering the presence of links. A first example is provided by the two 'friendship' networks where sparse subgraphs can be considered as communities whenever their connections are 'heavy' enough; a second example is provided by the weighted communities of Star Wars characters: the three major clusters, in this case, are the one induced by the characters of Episodes I-III, the one induced by the characters of Episodes IV-IX and the one contaning the villains of Episodes VII-IX.
the SBM, hence employing the information on the link density, both in a global and in a block-wise fashion.
Surprise becomes the best performer when binary, directed configurations are considered -see the second row of fig. 6: while the performance of modularity starts decreasing as soon as the value of µ t is risen and NMI Inf omap 0 when µ t crosses the value of 0.6, the performance of surprise 'degrades' much more slowlyin fact, for some instances of the LFR benchmark, it achieves a large value of NMI even for values µ t ≥ 0.8.
Let us now comment on the performance of our algorithms when weighted configurations are considered. The LFR benchmark, in fact, can be extended to account for weights as well. In this case, there are two mixing parameters: while the first one is the usual, topological mixing parameter, the second one, that will be indicated with µ w , accounts for the heterogeneity of weights -more precisely, controlling for the percentage of a node strength to be distributed on the links that connect it to the nodes outside its own community. The third parameter to be determined is the exponent of the strength distribution -that has been set to 1.5 for all realizations considered here. The results are, again, shown in fig. 6 (to be noticed that we have kept one of the two parameters fixed and studied the dependence of NMI and ARI on the other: specifically, we have frozen the topological mixing parameter and studied the dependence of the results Link weights refine the picture provided by just considering the presence of links: the core (in black) is, in fact, constituted by the nodes connected by the 'heavier' links, irrespectively from the link density of the former one. Bottom panels: results of the application of our framework for the detection of weighted 'bimodular' structures on the electronic Italian Interbank Money Market (e-MID) [42], for the maintenance period approximately corresponding to May 2009 (third panel: binary version; fourth panel: weighted version). As the picture shows, the (purely binary) core links are also the 'heavier' ones -an evidence confirming an ubiquitous tendency in economic and financial systems, i.e. binary and weighted quantities are closely related.
on µ w , thus inspecting the performance of our algorithms as the weights are redistributed on a fixed topology). Infomap is, again, a strong performer although its performance keeps decreasing abruptly as µ w exceeds a threshold value depending on the particular setting of the LFR benchmark; modularity, instead, performs worse than in the binary case although it is still more robust than Infomap. Although 'degrading' less sharply than Infomap, the performance of the purely weighted surprise seems to be the worst, here; on the other hand, the enhanced surprise outperforms the competing algorithms for intermediate values of the topological mixing parameter, irrespectively from the size of the communities: in fact, NMI surprise = 1 even for µ w = 0.8. Similar considerations hold true when weighted, directed configurations are considered (with the only difference that, now, modularity steadily performs worse than the other algorithms, except for the largest values of the mixing parameter). As for the binary cases, surprise competes with modularity although it employes much less information than the latter: in fact, while the benchmark employed by modularity now coincides with the (sparse version of the) Weighted Configuration Modelhence, encodes the information on the entire strength sequence -surprise compares the WRGM with the WSBM, hence employing the information on the magnitude of the total weight, both in a global and in a block-wise fashion.
Let us now consider the RC benchmark. It consists of 512 nodes, grouped in 16 communities arranged in a ring-like fashion and whose sizes obey a power-law whose exponent has been set to 1.8; the smallest community is composed by 3 nodes. Such a configuration is progressively 'degraded' according to the following mechanism: first, a percentage p of links is selected randomly and removed; afterwards, a percentage p of links is selected randomly and rewired. In other words, a single 'degradation' parameter p drives the evolution of the initial ringof-cliques towards a progressively less-defined clustered configuration.
Results on specific implementations of the RC benchmark are shown in fig. 7: surprise outperforms both competing algorithms across the entire domain of the 'degradation' parameter p. More specifically, while modularity 'degrades' slowly as the value of p is risen, Infomap 'degrades' abruptly as p ≥ 0.4. Hence, for small values of such a parameter, Infomap outperforms modularity; on the other hand, for large values of p, modularity outperforms Infomap (although both NMI modularity and ARI modularity achieve a value which is around 0.6, i.e. already far from the maximum). Interestingly, for small values of p, the performance of Infomap and that of surprise overlap, both achieving NMI and ARI values which are very close to 1: as p crosses the value of 0.4, however, the two trends become increasingly different with Infomap being outperformed by modularity which is, in turn, outperformed by surprise. From a more general perspective, these results confirm what has been already observed elsewhere [39], i.e. that the best-performing algorithms on the LFR benchmarks often perform poorly on the RC benchmarks and vice versa.
Let us now inspect the performance of surprise in recovering binary 'bimodular' structures. To this aim, we have defined a novel benchmark mimicking the philosophy of the RC one, i.e. progressively 'degrading' an initial, well-defined configuration: • let us consider N c core nodes and N p periphery nodes. The core is completely connected (i.e. the link density of the N c × N c block is 1) and the periphery is empty (i.e. the link density of the N p × N p block is 0). So far, our benchmark is reminiscent of a core-periphery structureà la Borgatti-Everett; • let us now focus on the topology of the N c × N p  [43]. Top panels: upon running W in a hierarchical fashion, a 'core-within-the-core' is detected, revealing that the richest world countries (i.e. Canada, USA, the richest European countries, China and Russia -right panel, in dark green) constitute an even more tightly-linked cluster of nations among the ones with the largest trading activity (left panel, in black). Bottom panels: E penalizes the countries with both a small degree and a small strength; in fact, the second run exlcudes the Russia, a result seemingly indicating that while its strength is 'large enough' to be a member of the core, its degree is not.
bipartite network embodying the connections between the core and the periphery: in particular, let us consider each entry of such an adjacency matrix and pose a cp = 1 with probability p cp . Upon doing so, such a subgraph will have a link density amounting precisely at p cp ; • let us now 'degrade' such an initial configuration, by progressively filling the periphery and emptying the core. This can be achieved by 1) considering all peripherical node pairs and link them with probability q; 2) considering all core node pairs and keep them linked with probability 1 − q (or, equivalently, disconnect them with probability q). Upon doing so, we end up with a core whose link density is precisely 1 − q and with a periphery whose link density is precisely q. Now, varying q in the interval [0, p cp ] allows us to span a range of configurations starting with Borgatti-Everett and ending with Erdös-Rényi.
Specifically, here we have considered N c = 100, N p = 300 and p cp = 0.5. The result of our exercise is shown in fig. 8: as expected, the performance of the surprise worsens as the 'degradation' parameter becomes closer to p cp = 0.5; however, both the NMI and the ARI indices steadily remain very close to 1 -a result meaning that surprise optimization is not only able to correctly classify true positives (i.e. to keep the nodes originally in the same communities together) but also the other, possible kinds of node pairs. As with the exercise on community detection, let us now ask us if the presence of weights affects the detection of mesoscale structures. Generally speaking, the answer is, again, yes. Let us consider a toy core-periphery network: rising the weight of any two links connecting the core with the periphery allows the two nodes originally part of the periphery to be detected as belonging to the core. Analogously, if a bipartite topology is modified by adding weights between some of the nodes belonging to the same layer, W will detect a core-periphery structure as significant, the core nodes being the ones linked by the 'heaviest' connections.
Testing surprise on real-world networks. Let us now apply our formalism to the detection of mesoscale structures in real-world networks.
When coming to study real systems, particularly insightful examples are provided by social networks. To this aim, let us consider the one induced by the co-occurrences of characters within the Star Wars saga (i.e. the three trilogies) [40]. As shown in fig. 9 we have both considered the binary and the weighted version of it. In both cases two major clusters are visible. For what concerns the binary version of such a network, the optimization of S reveals the presence of two major clusters: remarkably, those clusters are induced by the characters of Episodes I-III (e.g. Yoda, Qui-Gon, Obi-Wan, Anakin,  fig. 9). As the optimization of W reveals, while fully connected subsets of nodes are considered as communities in case links have unitary weights, sparser subgraphs can be considered as communities as well whenever their inner connections are 'heavy' enough. On the other hand, both bottom panels of fig. 9 seem to confirm that one of the main limitations of surprise-like functionals is that of recovering a large number of small cluster of nodes.
Let us now compare the performance of S and W in order to see if, and how, the presence of weights affects the 'bimodular' mesoscale organization of networks. To this aim, let us focus on the network of co-occurrences of the characters of the novel 'Les Miserables' [41]. As shown in fig. 10, link weights indeed modify the picture provided by just considering the simple presence of links (see also [12]): the core of the weighted network is, in fact, constituted by the nodes connected by the 'heavier' links, irrespectively from the link density of the former one.
After having applied our framework to the analysis of social networks, let us move to consider financial networks. One of the most popular examples of the kind is provided by the electronic Italian Interbank Money Market (e-MID) [42], depicted in fig. 10. Notice that, for such a network, the vast majority of core links are also the 'heavier' ones, an evidence confirming a tendency that is ubiquitous in financial and economic systems, i.e. binary and weighted quantities -even at the mesoscaleare closely related.
Remarkably, the surprise-based formalism presented in this paper can be employed in a hierarchical fashion to highlight either nested communities or nested 'bimodular' structures. To clarify this point, let us consider the World Trade Web (WTW) in the year 2000 as a casestudy [43]. First, let us run W to highlight the core portion of the weighted version of such a network; as fig.  11 shows, the bipartition distinguishes countries with a large strength from those whose trade volume is low (basically, a bunch of African, Asian and South-American countries). Repeating our analysis within the core portion of the network allows us to discover the presence of a (statistically-significant) nested core: in fact, the secound-round optimization reveals that the 'core-insidethe-core' is composed by countries such as Canada, USA, the richest European countries, China and Russia.
Let us now compare it with E , run in a hierarchical fashion as well. The results of our exercise are shown in fig. 11. As evident from looking at it, the enhanced surprise is more restrictive than the purely weighted one, as a consequence of constraining the degrees beside the strengths. Hence, while the first run excludes the countries with both a small degree and a small strength, the second run exlcudes the Russia, a result seemingly indicating that while its strength is 'large enough' to allow it to be a member of the core, its degree is not.
In a sense, the optimization of E corrects the picture provided by the optimization of W as the core becomes less populated by 'low degree' nodesan effect which is likely to become more evident on systems that are neither financial nor economic in nature.
Let us now run, and compare, modularity, Infomap and surprise on the bunch of real-world networks above. Table II sums up the results. A first observation concerns the number of detected communities: while Infomap is the algorithm producing the smallest number of clusters, surprise is the one producing the largest number of clusters -more precisely, surprise outputs more and smaller clusters than the other two methods. As a consequence, our three algorithms produce partitions with an overall small overlap, as indicated by the NMI; the ARI confirms such an observation -although indicating that the pictures provided by Infomap and surprise are (overall) more similar than those provided by modularity and surprise (and, as a consequence, by modularity and Infomap). Interestingly enough, the values of the AWI are quite large -and larger than the corresponding NMI and ARI values: since it just focuses on the percentage of true positives, a good performance under such an index indicates that the two tested algorithms 'agree' on the nodes to be clustered together (although they may not -and, in general, will not -agree on the number of communities). Hence, the discrepancy between the ARI and the AWI may be explained by the presence of statistical 'noise' (i.e. 'misclassified' pairs of nodes, although the word may not be correct as the information about the 'true' partition is not available) around the bulk of nodes to be put together.
Let us stress once more that, whenever real-world networks are considered, information about the existence of a 'true' partition is rarely available; for this reason, exercises as the one we have carried out here may be useful to gain insight on the system under study: instead of trusting just one algorithm, combining pairs of theme.g. by considering as communities the subsets of nodes output by both -may be the right solution to overcome the limitations affecting each single method.
The 'SurpriseMeMore' Python package. As an additional result, we release a comprehensive package, coded in Python, that implements all aforementioned surprise-like variants.

CONCLUSIONS
The hypergeometric distribution -together with its many variants -has recently revived the interest of researchers who have employed it to define novel network ensembles [45], recipes for projecting bipartite networks [46], etc.
Remarkably, the distributions related to it allow for a wide variety of benchmarks to be defined, each one embodying a different set of constraints. In the present paper we have explored the 'power' of the hypergeometric formalism to carry out the detection of mesoscale structures: remarkably, it allows proper statistical tests to be definable for revealing the presence of modular, coreperiphery and bipartite structures on any kind of network, be it binary or weighted, undirected or directed. According to the classification proposed in the introductory section of the paper, we believe surprise to belong to the second class of algorithms -its asymptotic expression embodying a sort of LRT aimed at choosing between alternative hypotheses (the RGM and the SBM, the WRG and the WSBM, etc.).
More in general, our approach reveals the superiority of the algorithms for mesoscale structures detection belonging to the second and to the third class with respect to those belonging to the first one: still, the two classes of statistically-grounded approaches compete on some, specific benchmarks, as the comparisons carried out on the LFR and the RC ones clearly show (specifically, methods performing well on LFR benchmarks do not on RC benchmarks and vice versa). This also suggests a strategy to handle mesoscale structures detection on real-world networks: as the information on the presence of a possible, 'true' partition is rarely available, a good strategy may be that of running different algorithms on the same empirical networks, check the consistency of their output and combine them, e.g. by taking the overlap -in a way that is reminiscent of multi-model inference.
Although the surprise-based approach is powerful and versatile, its downside is represented by the specific kinds of tests that are induced by the optimization of surpriselike score functions: comparisons between benchmarks 'ignoring' the local structure of nodes (i.e. their degree) are, in fact, carried out. While this seems to be perfectly reasonable when considering core-periphery structuressee also the contribution [29] whose authors claim that a core-periphery structure is always compatible with a network degree sequence -this is no longer true for the community detection task [26] -indeed, as it has been noticed elsewhere, ignoring degrees may be at the origin of the large number of singletons output by surprise as assigning nodes with few neighbors to larger clusters may be disfavored from a statistical point of view.
The observations above call for the 'extension' of the hypergeometric formalism we have studied here to include more refined benchmarks as the ones constraining the entire degree sequence. As stressed in the main text, community detection on binary networks can be implemented by carrying out an exact statistical test whose function f must be identified with a binomial hypergeometric distribution. Hence, we can compactly write with a clear meaning of the symbols. The correctedness of such the definition of can be explicitly checked upon calculating its normalization via the Vandermonde's identity The generalization of the formalism above to detect 'bimodular' structure on binary networks is straightforward. It is enough to consider the multivariate (or multinomial) version of the hypergeometric distribution, i.e.
the correctedness of such a definition of can be checked upon calculating its normalization by applying the identity The asymptotic expression of S and S can be derived by applying the Stirling approximation n! √ 2πn n e n to the binomial coefficients appearing into the definition of the versions of the surprise described above, leading to expressions like where Ber(x, y, z) = z y (1 − z) x−y is a Bernoulli probability mass function, z = y x and p = L V .
The time, in seconds, to find the corresponding optimal partitions is reported in tables III and IV: overall, denser networks require more time.  When analyzing weighted networks, one needs to move from the binomial hypergeometric distribution to its negative counterpart, i.e. to with a clear meaning of the symbols. The soundness of the definition above can be cheked upon calculating the normalization of the distribution NH(w • |V + W, W, V • ) via the identity The position above induces the definition of the weighted surprise as The generalization of the formalism above to detect 'bimodular' structure on weighted networks is straightforward. It is enough to consider the multivariate (or multinomial) version of the negative hypergeometric distribution, i.e.
whose correctedness can be checked upon calculating its normalization, by applying the identity The asymptotic expression of W and W can be derived by applying the Stirling approximation n! √ 2πn n e n to the binomial coefficients appearing into the definition of the versions of the surprise described above, leading to expressions like where Geo(x, y, z) = z y (1 − z) x is a geometric probability mass function, z = y x+y and q W V +W .
The time, in seconds, to find the corresponding optimal partitions is reported in tables III and IV: overall, denser networks require more time.  The enhanced hypergeometric distribution is a negative hypergeometric distribution 'conditional' to the existence of connections within communities, i.e. l • > 0. It reads the first expression being valid when l • < L and the second expression being valid in the 'extreme' case l • = L (i.e. whenever the total weight has to be redistributed on top of links which are all put inside the communities). Normalization is ensured by the fact that where we have used the relationships as it should be: in fact, the expected value of (the number of) connections should not depend on the weights put on the connections. Moreover, The asymptotic expression of E can be derived by applying the Stirling approximation n! √ 2πn n e n to the binomial coefficients appearing into the definition of the enhanced surprise, leading to expressions like Let us now briefly consider the distribution inducing the definition of the enhanced bimodular surprise, i.e. the 'multivariate' enhanced hypergeometric distribution. As before, it is 'conditional' to the existence of connections within either modules one is looking at (i.e. either l • > 0 or l • > 0) and reads with a clear meaning of the symbols. Notice that the cases l • = 0, 0 < l • < L and l • = 0, 0 < l • < L can be easily recovered from the first row of the definition above.

Appendix D: Pseudocodes for surprise optimization
Let us now provide a description of the algorithms implemented in the 'SurpriseMeMore' Python package (freely downloadable at the URL https://github.com/EmilianoMarchese/SurpriseMeMore) for surprise optimization: to this aim, we will also show the corresponding pseudocodes. Since an exhaustive exploration of the space of all possible partitions of a network is not feasible (especially when dealing with large graphs), one often proceeds heuristically, in either an agglomerative fashion or fixing the number of clusters to be detected a priori.
Community detection on binary and weighted networks. The algorithm for community detection on binary networks is called PACO (PArtitioning Cost Optimization) and its pseudocode can be found in [15]. The algorithm implemented there is agglomerative in nature: in the initial configuration of the algorithm, in fact, nodes are treated as singletons (hence, there are as many communities as nodes). Then, edges are sorted according to the value of their Jaccard index, in a decreasing order: as the Jaccard index of an edge is proportional to the number of common neighbours of the two connected nodes, a large value of it likely points out that the two nodes must be assigned to the same module. For further details about PACO, see [15].
Here we adopt a similar strategy to the one implemented by PACO in the binary case, with two main differences. The first one concerns the edge-sorting step: edges, in fact, are no longer sorted according to the Jaccard index but randomly ordered; this strategy avoids getting stuck in local minima more effectively. The second one concerns the possibility of letting entire communities merge (instead of moving single nodes from one to another) with a certain probability p mix . if 0 ≤ unif ormDistribution(0, 1) ≤ p mix then 17: for node t ∈ C[u] do  ⇒ consider each node and swap its membership with the one of each first neighbour of it; accept the move if S , W or E decreases; 32: end for 33: ⇒ repeat the for-loop to improve the chance of finding the optimal partition Remarkably, it is also possible to define a version of the surprise-based community-detection algorithm, where the number of clusters to be detected is fixed a priori:  ⇒ consider each node and swap its membership with the one of each first neighbour of it; accept the move if S , W or E decreases; 26: end for 27: ⇒ repeat the for-loop to improve the chance of finding the optimal partition 'Bimodular' structures detection on binary and weighted networks. In this case, the algorithm we implement fixes the number of clusters to be detected a priori, since each node of the network can belong to only one out of two possible subsets (e.g. the core and the periphery) -a circumstance that is reflected into the initialization of the algorithm, i.e. a vector whose entries are randomly chosen to be either 0 or 1. A second difference with respect to the community detection case concerns the criterion according to which edges are sorted (for purely bipartite graphs, in fact, the Jaccard index is zero for all edges, as nodes connected by an edge always lie on different layers): in our modified version of the PACO algorithm, when binary graphs are considered, we sort links according to their eigenvector centrality, which has been shown to proxy 'coreness' to quite a good extent. As a final step, we consider a number of random re-assignments of nodes with the aim of preventing the possibility of getting stuck in a local minimum (a random move consists of selecting 3 nodes belonging to the same group and evaluating if assigning them to the other group would further minimize surprise).
The edge-sorting strategy described above is not the best one in the weighted case: in this case, in fact, randomlyordered edges are to be preferred, as they avoid getting stuck in local minima more effectively.   ⇒ randomly switch the membership of 3 nodes in the same group and accept the move if S , W or E decreases; 23: end for 24: ⇒ repeat the for-loop to improve the chance of finding the optimal partition