Revealing the predictability of intrinsic structure in complex networks

Structure prediction is an important and widely studied problem in network science and machine learning, finding its applications in various fields. Despite the significant progress in prediction algorithms, the fundamental predictability of structures remains unclear, as networks’ complex underlying formation dynamics are usually unobserved or difficult to describe. As such, there has been a lack of theoretical guidance on the practical development of algorithms for their absolute performances. Here, for the first time, we find that the normalized shortest compression length of a network structure can directly assess the structure predictability. Specifically, shorter binary string length from compression leads to higher structure predictability. We also analytically derive the origin of this linear relationship in artificial random networks. In addition, our finding leads to analytical results quantifying maximum prediction accuracy, and allows the estimation of the network dataset potential values through the size of the compressed network data file.

P redicting structure or links of network is commonly defined as estimating the likelihood of existence of unobserved links or potential new links [1][2][3][4][5] . Network as a common form of data representation is ubiquitous across a broad range of fields from biology [6][7][8] , recommendation systems 9,10 to social media 1,3,11 . It represents the complex relationships or interactions among the elements in a system, which usually cannot be well described by simple mathematics. Hence, machine learning has been widely used in link predictions [10][11][12][13][14][15][16][17][18] . For instance, predicting protein-protein 7,8 or drug-target [19][20][21] interactions can guide more accurate biological experiments and reduce the experimental costs and time 2,6 . Despite the immense ongoing efforts in developing prediction algorithms, the fundamental understanding of the intrinsic prediction limit that can provide the much needed guidance is lacking. The difficulty lies with the fact that it is almost impossible to know the exact underlying mechanism of the network formation. In addition, the real networks are usually highly complex and gigantic in size, with many short feedback loops which cannot be effectively analyzed-a challenge 22 faced by both statistical physicists and computer scientists. Hence, understanding the prediction limits and quantifying it remains to be a long-standing challenge 2,23 .
In this study, we reveal the intrinsic predictability of networks through their structural properties. Intuitively, a network structure that can be captured with a few words means it is simple, and its links are easily predictable such as one-dimensional chain and two-dimensional lattice. Conversely, if a network requires lengthy description, then it has very complicated structure and its links are hard to predict. In computer language, the structure of any network can be encoded into binary strings. This motivates us to find the underlying relationship between the network compression length of the shortest binary string and the prediction limit. On the other hand, the inherent prediction limit is the maximum predictability or performance that the theoretical best predicting algorithm (TBPA) can achieve. Here the network structure predictability is defined as the performance of TBPA and is quantitatively measured using entropy in this paper. Thus, without knowing the exact underlying dynamics of network formation which determines this limit, we can use the best prediction algorithm available (BPAA) 1,2,23,24 to approximate the performance of TBPA. With these two quantities of the shortest lossless compression length and performance of TBPA, here, for the first time, we discover a linear relationship between them in different empirical networks such as biological networks, social networks, and technology networks (see Supplementary Note 1 for the detailed description of studied networks). Our finding implies that the shortest compression length of a network can tell us the structure predictability, which sets the limit of any prediction algorithm.

Results
Network shortest compression length. The shortest possible compression length can be calculated by a lossless compression algorithm 25 , which is a proven optimal compression for random networks and efficient for many real networks. Since all of the structural properties of two isomorphic networks are exactly the same, the algorithm first encodes the network structure into a string of binary codes, from which an isomorphic network can be reconstructed. To further remove the correlations in the structure, the binary string is compressed by a recursive arithmetic encoder 26 , which exploits the dependencies between the symbols in this binary string. After these two compression operations (see Fig. 1a), the length of the final bit string should be close to the network's structural entropy and can well measure its randomness (Shannon's source coding theorem 27 , see Supplementary Note 2). The length of this bit string is expected to increase as the structure becomes more random, and has been validated in Fig. 1b, which shows shuffled networks has longer compression length compared with the original empirical networks. As randomness in shuffling increases, the compression length increases monotonically.
Naturally the compression length is longer for a network with more nodes and links given the same randomness, i.e., the size of the network contributes to the compression length, rather than just the level of randomness of the network itself. In order to remove the size effect of the network, we normalize the compression length L through dividing it by the theoretical maximum compression length R, which corresponds to Erdős-Rényi (ER) 25,28 network of the same number of nodes and links (see Supplementary Note 2). The normalized value of the shortest compression length L * is given by Here R ¼ N 2 À Á hðqÞ À NlogN 25 , where N is the number of nodes and q is the probability having a link between any pair of nodes, also expressed as the binary entropy given by hðqÞ ¼ Àqlogq À ð1 À qÞlogð1 À qÞ.
Here log denotes the logarithmic operation with base 2, and we use this convention throughout this work.
Structure predictability measured by algorithm performance. Now we compare the shortest compression length with the performance of BPAA which approximates the structure predictability of network. In a large literature on link prediction algorithm [1][2][3][4][5][6]11,15,16 , it is common to assign each unlinked pair of nodes (possible missing link) a score, and higher score means higher likelihood of this pair being a missing link. Here we adopt the leave-one-out 31,32 approach and use the score of the existing links to quantify the BPAA performance. First we remove a link e i , and then use one particular prediction algorithm to estimate a score for each of the unlinked node pairs including e i . Based on the ranking of the scores in descending order, we get the ranking r i of the removed link e i . So when r i = 1, it means that the algorithm tells us the removed link is the most probable missing link among all other unlinked node pairs. We carry out this calculation of r i from the original network for every link at a time, and obtain a sequence of rank positions D = {r 1 , r 2 , . . . , r E }, where E is the total number of links in the original network. Naturally, the entropy of the distribution of D is a good holistic measure of the algorithm performance. For example, an ideal algorithm for a highly predictable network would have D = {r 1 = r 2 = ⋯ = r E = 1}, thus yielding the lowest distribution entropy of r i . Conversely, an ideal algorithm for a network with low predictability has very different values in D, leading to a high entropy of its distribution. In our calculation of this algorithm performance entropy H, the value r i can vary in the range 1 r i NðNÀ1Þ 2 À hkiN 2 þ 1 % N 2 2 , where NðNÀ1Þ 2 À hkiN 2 þ 1 is the total number of unlinked pairs, and 〈k〉 is the network's average degree. Thus we divide such range into bins with equal width N to avoid the contribution of network size N on the result (see Supplementary Note 3 for the discussion of other bin widths), and calculate H based on the probability distributions of the N∕2 bins: H ¼ À P N=2 j¼1 p j logp j where p j is the probability of r i in bin j (Supplementary Note 3). Figure 1c illustrates an example of such distribution for the Metabolic network 29 by RA algorithm 30 .
It is worth mentioning that using rank distribution entropy to measure the algorithm performance is the innovation in our work. In particular, we implement the leave-one-out method 31,32 to obtain the rank distribution. The other ranking methods 2,3,23 which pre-remove x fraction of links (x is usually 5% or 10%) with multiple samples are similar to the leave-one-out method, and would yield similar results as leave-one-out if x is very close to 0 (see Supplementary Note 4 and Supplementary Fig. 6). Using the leave-one-out method, there are three main advantages. Firstly, it is parameter-free, i.e., we do not need to consider the impact of the choice of x, and allows analytical study. Secondly, doing so preserves the original network structure as much as possible, such that the intrinsic true predictability is also preserved. Thirdly, the leave-one-out method has negligibly small fluctuation (see Supplementary Fig. 8). The main drawback is the higher computational complexity compared with only using 1 sample of removing x fraction of links. However, if multiple samples of removing x fractions are considered, the advantage of complexity of above other methods will not be significant over ours.
In order to obtain a good approximation of the TBPA performance of the network that is less dependent on a particular prediction algorithm being used, we employ a range of widely applied 11 prediction algorithms (see Supplementary Note 5), such as structural perturbation method (SPM) 23 , local random walk (LRW) 33 , average commute time 24 and common neighbor 1 to calculate the ranking positions r i , and use the one that gives the lowest H value H BPAA (Table 1) as the closest estimate of the network's predictability H TBPA . To further explore the relationship between network structure and H BPAA , it can be seen in Fig. 1d that as the network structure is progressively randomized due to shuffling, the entropy H BPAA increases, signifying a decrease in predictability. Eventually with enough shuffling, the algorithm performance entropy H approaches that of ER networks with a nearly uniform distribution of D (Fig. 1d inset). To further remove the effect of network size N and average degree 〈k〉, we normalize the BPAA performance entropy H BPAA value by logN À 1, which is the BPAA entropy of an ER network of the same size and average degree. Hence the normalized BPAA performance entropy H BPAA is defined as Empirical linear relationship. Having both network's structural shortest compression length and the BPAA performance, we are able to find the relationship between the two. Surprisingly, we discover a clear linear relationship between the H Ã BPAA and L * across 12 empirical networks in very different fields including biology, social media, transport and economics as shown in Fig. 2a: It is worth mentioning that, the way H and L are normalized is critical to obtain the linear relationship, as without proper normalization, the network's compression length cannot serve as a good indicator to quantify the structure predictability of a network (see Supplementary Note 6).
We find that such a linear relationship in Eq. (3) still exists even after we shuffle links as seen in Fig. 2b, or randomly add/ delete links as seen in Fig. 2c property of the shortest compression length 25,27 , i.e., the sum of the shortest compression lengths L A , L B of two independent networks A and B (with the same nodes but different links) is the same as the compression length L A+B of the combined network A + B (Fig. 2e). This is also valid for the compression algorithm in case of shuffling and adding links at random on empirical networks ( Fig. 2f, g). We now illustrate how to use the linear relationship in Eq. (3) to quantify the performance of a prediction algorithm. By comparing the Jaccard algorithm 35 performance entropy with the linear relationship (see Supplementary Fig. 11), it can be seen that for a given network, the further is the performance entropy from the linear line, the less optimal it is from the best possible prediction results (see Supplementary Note 7). This quantifiable distance can serve as a benchmark to practical algorithms in a given network. In other words, the further it is from the linear line, the more potential improvement it is possible for any new algorithm to achieve, warranting further effort on algorithm development. Conversely, if the algorithm performance entropy lies close to the line, it means that the algorithm is already performing quite well, and further exploration on algorithms is unlikely to yield significant improvement, suggesting diminishing returns in developing new algorithms.
Since the compression process can be executed through different node sequences, we also study the effect of compressed node sequences on the compression length, such as random, high degree priority and low degree priority. We find that the differences in compression lengths are very small (see Supplementary Fig. 2). Hence, we can see this compression algorithm is reliable when operating on real networks. Naturally, such relationship between H Ã BPAA and L * indicates that the link prediction limit of a network can be directly inferred from the network topological structures. More explicitly, by compressing the network structure into binary string and calculate the string's length, one can estimate the link prediction limit H Ã BPAA . Apart from the theoretical significance of the relation, one practical advantage of estimating H Ã BPAA with L * is that L * is only dependent on the network topological structure and has low computational complexity O(N + E) 25 . In contrast, estimation of H Ã BPAA requires a large number of available prediction algorithms, which usually have high computational complexity. Specifically, for calculating the rank distribution entropy H, we need to employ an algorithm for every link in a network to obtain its ranking. Consequently, the computational complexities of algorithms like SPM 23 and LRW 33 are about O(N 3 E) and O(N〈k〉 n E) in estimating H, where n is an integer constant. Furthermore, an additional benefit of the compression length is that it can be served as an independent indicator to identify missing or false links (see Supplementary Note 8).
Theoretical linear relationship. The empirical linear relationship inspires us to further figure out the underlying mathematical connection between the network shortest compression length and link prediction limit. For simplicity, we assume that an artificial network is generated from a static random matrix Q whose entry q ij denotes the link formation probability between node i and j. According to Shannon's source coding theorem 27 , the shortest compression length L of this artificial network is L ¼ P i>j hðq ij Þ À NlogN, which is the structural entropy 25 , where h(q ij ) is the binary entropy given by hðq ij Þ ¼ Àq ij logq ij À ð1 À q ij Þlogð1 À q ij Þ. To simplify representation, we use U to denote the term À P i > j q ij logq ij and expand the logarithm as a Taylor series, of which the higher-order small terms can be neglected, yielding: where ln denotes the natural logarithm, and we use this convention throughout this work. Since the occurrence of each link in this artificial network solely depends on the q ij , it is certainly true that the most accurate score that a TBPA could achieve should be exactly equal or proportional to q ij . Thus, without employing any prediction algorithm, we can obtain the ranking sequence D (see the right part of Fig. 3a) directly from Q and are able to numerically quantify H Ã TBPA in the same way that we calculate H Ã BPAA . We find that H Ã TBPA can be theoretically estimated by Q, given that the TBPA ranking distribution clearly agrees well with the distribution of q ij in Q (Fig. 3b). The entropy of the latter is given by Recalling that in the calculation of H BPAA , we have divided the ranking range into bins with equal width N. In such case, this coarse-grained distribution of Q is done by replacing every N values of q ij (in the descending order) by their average valueq ij (see the left part of Fig. 3a), and the  resulting entropy satisfiesH Q ¼ H Q À logN (see Supplementary Note 9). The difference between H TBPA andH Q is due to the fact that for H TBPA , link prediction algorithms only measure the likelihood of unobserved pairs of nodes forming a link. Such difference, however, contributes negligibly to the calculation of entropy as shown in Fig. 3b, yielding H TBPA %H Q . Taken together, we obtain (see details in Supplementary Note 9): We then normalize H TBPA by logN À 1 through the same way of Eq. (2), yielding Combining with Eqns. (1), (4)-(6) and eliminating the variable U, we obtain the linear relationship between L * and H Ã TBPA : When N is large Eq. (7) simplifies to In the thermodynamics limit Eq. (8) can be further approximated when logN ) loghki: Usually for many real networks, loghki logN is not negligible, and Eq. (8) is a better approximation. The detailed mathematics above is given in Supplementary Note 9.
To validate Eq. (7), we first construct artificial networks based on each empirical network in Fig. 2 with the same degree sequence, but a simple network formation mechanism based that Matrix Q distribution 5 6 TPBA performance entropy, H TBPA 7 8 9 Fig. 3 Theoretical relationship between the shortest compression length and network prediction limit. Networks are generated from a static random matrix Q with elements generated from the degree distribution of empirical networks. a Comparison between the ranking distribution from TBPA and that from random matrix Q based on Metabolic network 29 . The left half represents the ranking distribution of coarse-grained probabilitiesq ij , each of which is the average value of every N values of q ij ranked in a descending order. The right half illustrates the TBPA ranking distribution of existing links for this artificially generated network. The probability p j is an average value from 100 simulations. b Values of coarse-grained entropyH Q from the probability distributionq ij vs. TBPA performance entropy H TBPA for all of the artificial networks generated from the empirical ones. Each dot corresponds to the value of (H TBPA ;H Q ) of an artificial network, showing a good match between the two values. c Theoretical linear relationship between L * and H Ã TBPA calculated based on the empirical networks' degree distributions. Points of the same shape and color correspond to an empirical network (bottom-left) and its different shuffled versions, i.e., having its links shuffled different number of times up to all of original links. Each point represents the value pair of ðL Ã ; H Ã TBPA Þ of an artificial network with the same degree distribution as a corresponding (shuffled) empirical network. L * is calculated from Eq. (1) and H Ã TBPA is calculated from Eq. (6). The blue solid line is the average of 12 studied networks' analytical slopes obtained from Eq. (7) and the shaded region denotes the standard deviation. The gray dash line is the analytical result of Eq. (9) for the limit of 〈k〉 → ∞. d Competition between two terms of the slope in Eq. (8). e The empirical and theoretical value of slope in the linear relationship between L * and H Ã TBPA . The purple plane denotes the plane of empirical value at 1.63, and each colored point on the plane represents a real network. The lower curved surface below represents the theoretical values of the slope given by Eq. (7), and each colored point on the lower surface represents an artificial network constructed with the same degree distribution as original empirical networks, and the subsequent points with the same color are the artificial networks with randomly added links to increase 〈k〉. The red curve (lnN ¼ hki) on the lower surface is an estimated boundary (see Supplementary Note 9) that our theories are not valid far from the left side of it.
only depends on the static probability distribution matrix Q, with each element q ij ¼ k i k j 2E , where k i denotes the degree of node i in an empirical network. Secondly, we shuffle each above network by randomly picking a fraction of links and rewiring them randomly (same operation in Fig. 2b) and then construct artificial networks based on each shuffled empirical network's degree sequence. In the shuffling process, the randomness of the network increases with the increasing of the shuffling strength. Theoretically, the values (H Ã TBPA , L * ) changes along the straight line predicted by Eq. (7). It can be seen from Fig. 3c that this theoretical relation in Eq. (7) is well captured by our simulation for the different networks and their shuffled versions. We have also considered two other genuine edge-independent synthetic networks: the degree-correlated stochastic block model 36 characterizing the community structure, and the latent-geometric network model 37 characterizing network spatial structure. Both models have been shown to reproduce certain structural properties of real network. We find that the structure predictability properties obtained by these two models also follow our analytical result for artificial networks (see Supplementary Note 9 and Supplementary Fig. 15).
Surprisingly, the artificial networks generated by matrix Q leads to ðL Ã ; H Ã BPAA Þ pairs falling on a straight line as seen in Fig. 3c, with a slope that is different from the one given by Eq. (9) from the thermodynamics limit approximation. A careful examination shows that the artificial networks' average degrees 〈k〉 increase with the network sizes N faster than logN, balancing out the changes in both loghki logN and 2 hki in Eq. (8) (see Fig. 3d), therefore keeping the slope relatively constant around 0.7 instead of the thermodynamics limit 1 (see the plateau in Fig. 3e). For the approximation result in Eq. (9), we find that the thermodynamic limit slope 1 with large 〈k〉 can be realized in the artificial network with size larger than 10 100 (see Supplementary Fig. 14). But this large number is almost impossible in real networks. Moreover, the values of the slope between H Ã TBPA and L * of artificial networks are significantly lower than the empirical slope of 1.63 (Fig. 3e). This may be due to the more complex mechanisms and constraints in empirical network formations compared with simplistic random link connection in the artificial networks. Intuitively one expects the structure predictability of the empirical network to be higher than that of the purely random due to such additional constraints or mechanisms, and it is reflected in higher slope observed. Such hypothesis does not directly address the differences quantitatively, yet it sheds some light on the complex relationship between network entropy and structure predictability. At the same time, the empirical value of 1.63 hints at some strong universal mechanism that is common to all of these empirical networks.
Bounds of link prediction precision. In many practical applications involving link prediction, like recommendation system, a prediction algorithm gives the most likely missing links in the network among all possible links. One way to assess the algorithm's prediction precision is through the success of prediction [1][2][3][4][5][6]11,15,16 . In other words, we can look at the probability p 1 that the distribution of ranks given by the algorithm on the removed links defined in D. Note that if we take the top N predicted links to calculate the standard precision, then it is equivalent to p 1 used in this paper, since p 1 refers to the fraction of correct links among the top N predicted links. In practice, a good algorithm will give a ranking distribution that is usually monotonically decreasing, i.e., p 1 ≥ p 2 ≥ ⋯ ≥ p N∕2 ≥ 0 similar to Fig. 1b. That means the link predicted by the algorithm is indeed the most likely to be missing. Actually, we observe this phenomenon for all the algorithms and data used in this paper (see Supplementary Fig. 26). Using the linear relationship between L * and H Ã BPAA , for any network, we can calculate its normalized shortest compression length L * from its compressed binary string. Together with the assumption p 1 ≥ p 2 ≥ ⋯ ≥ p N∕2 ≥ 0, we arrive at an implicit function of the precision upper bound p 1 : The exact value of p 1 can be obtained by solving Eq. (10) (see Fig. 4a). Additionally, the lower bound p 1 of the prediction precision p 1 can be simply calculated through the following explicit formula (details in Methods section): One can define the prediction precision in a more general way other than p 1 . For instance, rather than defining the probability p 1 of the removed link falling into the first interval, one can loosen the definition to the first C intervals (1 ≤ C ≤ N∕2), i.e., defining precision as P C ¼ P C j¼1 p j (see Supplementary Note 10 for details on the derivations). Indeed, we validate the above upper bounds in the empirical networks as shown in Fig. 4b Commercial value of network dataset. In practice, the prediction bounds enable something that was not possible before, for instance it can be used to estimate the commercial value of a network dataset, through its compressed size without developing any prediction algorithm. Thus, for a compressed network with L bit, the conservative commercial value V is where Θ is an external economic variable. For example, in the scenario of inferring interactions between proteins, Θ can be written as θn, where θ represents the unit experimental costs that can be saved if successfully predicting one interaction and n(≤N) denotes the number of predicted interactions. In other words, once we obtain the smallest number of bits L needed to store a compressed network data file, we can directly derive the approximate value of this dataset through Eq. (12). Note that the commercial value we refer to above is the potential additional value that can be realized by predicting unobserved information from the network, without considering the overlapping value from external information outside the network structure. A more general framework of data commercial values when various external information is available is discussed in Supplementary Note 11.

Discussion
Although we have seen that our finding applies to a broad range of empirical networks and their artificial counterparts, it is important to take note of certain limitations in practice. It is know from ref. 38 that Eq. (4) is accurate when lnN ( hki ( N À lnN (see the red curve lnN ¼ hki in Fig. 3d). Therefore when a network is extremely sparse or dense, the relationship between L * and H Ã BPAA may not be close to what we found. In addition, the network entropy estimation and structure predictability analysis assume randomness in structure. That means for regular networks like lattice structure, our finding does not hold. To illustrate the impact of regular structural features, we combine a regular network from a circular model network 39 (Fig. 5a) and real networks. The analysis on such synthetic networks show that indeed, such regular structure in a random network makes the H Ã BPAA vs. L * relation deviates from the slope of 1.63 as seen in Fig. 5b. Hence our result is more valid for networks without significant amount of regular links. That being said, if a network is very regular in structure, lattice for instance, its structure is easy to predicted and compressed due to high regularity involved.
In conclusion, for the first time we established a theoretical framework to quantify the intrinsic structure predictability in real and artificial random networks solely based on network structure, and independent from any prediction algorithms. With theoretical intuition on proper normalization, we uncover a universal linear relationship between the shortest compression length of the network structure and its structure predictability for a wide range of real networks such as biological networks, social networks and infrastructure networks, and analytically derived the mathematical origin of this linear relationship in artificial networks. In principle, such relationship can serve as a benchmark to quantify the performance of any practical prediction algorithm in nonregular networks. Leveraging upon this linear relationship, we can obtain the structure predictability that is intrinsic to the structure of complex networks and provide the accuracy bounds for any link prediction algorithm. In practice, our method can also be used to estimate the commercial value of a network dataset through its compressed length without using any link prediction algorithm. Our finding is demonstrated upon data structure in the form of networks. However, if the structural entropy and prediction limit are linked through the underlying dynamical process for data structures other than networks, it is likely that the predictability in some other machine learning problems for different types of systems can also be inferred through similar approaches of optimal compression.

Methods
Compression algorithm. Our compression scheme builds upon the seminal paper of 25 , which is a two-step lossless compression of graphs. First we encode a network into two binary sequences B 1 , B 2 , through traversing all nodes on the network from an initial node, and encode nodes' neighbors according to specific rules 25 . In the second stage, both B 1 , B 2 are compressed by an improved arithmetic encoder based on recursive splitting proposed by K. Skretting 26 to exploit the dependencies between the symbols of the sequence. After that we obtain two compressed binary sequencesB 1 ,B 2 , and we define the compression length of the network L to be the total length of these two sequences: where 'ðB 1 Þ,'ðB 2 Þ are the length ofB 1 andB 2 , respectively. More details about the algorithm is provided in Supplementary Note 2.
Upper and Lower Bounds of p 1 and P C . Here we demonstrate how to identify the upper and lower bounds of link prediction precision p 1 and P C . The basic idea is to transform these problems into optimization problems of certain boundary conditions. Firstly, we calculate the network's BPAA performance entropy H Ã BPAA through its shortest compression length L * given by Eq. (3). Its un-normalized value is ð1:63L Ã À 0:63ÞðlogN À 1Þ, and by the definition of entropy it can be written as There are two constraints for p i s. The first is that they sum up to unity, i.e.   and: as they are arranged in decreasing order. With these boundary conditions, p 1 is maximized when the other probabilities are equal. Then it can be found that the upper bound of p 1 is given by Eq. (10).
And the minimum value of p 1 is when there are as many p i s as close to p 1 as possible, leading to a lower bound given by Eq. (11).
For predictability of the top C intervals P C ¼ P C j¼1 p j , the upper bound corresponds to the case when the last N∕2 − C probabilities are the same, yielding ÀP C log P C C À ð1 À P C Þ log 1 À P C N=2 À C ¼ ð1:63L Ã À 0:63ÞðlogN À 1Þ: For the lower bound of P C , the value of p C can be introduced to construct the minimization problem. It has the boundary condition of 0 ≤ p C ≤ P C ∕ C. The approximate lower bound is the solution of the following: P C % minfP C : ÀγðP C ; p C ÞlogγðP C ; p C Þ À ð1 À γðP C ; p C ÞÞlog p C ¼ ð1:63L Ã À 0:63ÞðlogN À 1Þg: The more detailed mathematics is provided in Supplementary Note 10.

Data availability
Data of this study is available at http://www.huyanqing.com/.